Diagnosis distribution analysis

library(here)
library(cowplot)
source(here("utils/data_processing.R"))
source(here("utils/figures.R"))

Import data

# Original outputs
df_gpt3.5 <- read_csv(here("data/processed_diagnoses/diagnoses_gpt3.csv.gz"))
df_gpt4.0 <- read_csv(here("data/processed_diagnoses/diagnoses_gpt4.csv.gz"))
df_claude3_haiku_t1.0 <- read_csv(here("data/processed_diagnoses/diagnoses_claude3_haiku_t1.0.csv.gz"))
df_claude3_haiku_t0.1 <- read_csv(here("data/processed_diagnoses/diagnoses_claude3_haiku_t0.1.csv.gz"))
df_claude3_opus_t1.0 <- read_csv(here("data/processed_diagnoses/diagnoses_claude3_opus_t1.0.csv.gz"))
df_gemini1.0_pro_t1.0 <- read_csv(here("data/processed_diagnoses/diagnoses_gemini1.0_pro_t1.0.csv.gz"))
df_gemini1.5_flash_t1.0 <- read_csv(here("data/processed_diagnoses/diagnoses_gemini1.0_pro_t1.0.csv.gz"))

# ICD mapped outputs
#TODO Find source of NAs 
#TODO drop index column in processing script
df_gpt3.5_icd <- read_csv(here("data/processed_diagnoses/diagnoses_gpt3.5_icd.csv.gz")) %>% drop_na() %>% select(-index)
df_gpt4.0_icd <- read_csv(here("data/processed_diagnoses/diagnoses_gpt4_icd.csv.gz")) %>% drop_na() %>% select(-index)
df_claude3_haiku_t1.0_icd <- read_csv(here("data/processed_diagnoses/diagnoses_claude3_haiku_t1.0_icd.csv.gz")) %>% drop_na() %>% select(-index)
df_claude3_opus_t1.0_icd <- read_csv(here("data/processed_diagnoses/diagnoses_claude3_opus_t1.0_icd.csv.gz")) %>% drop_na() %>% select(-index)
df_gemini1.0_pro_t1.0_icd <- read_csv(here("data/processed_diagnoses/diagnoses_gemini1.0_pro_t1.0_icd.csv.gz")) %>% drop_na() %>% select(-index)

Rank abundance

Individual model data

Original responses

rank_abundance_plot(df_gpt3.5)+ggtitle("ChatGPT 3.5")

rank_abundance_plot(df_gpt4.0)+ggtitle("ChatGPT 4.0")

rank_abundance_plot(df_claude3_haiku_t1.0)+ggtitle("Claude3 Haiku t1.0")

rank_abundance_plot(df_claude3_haiku_t0.1)+ggtitle("Claude3 Haiku t0.1")

rank_abundance_plot(df_claude3_opus_t1.0)+ggtitle("Claude3 Opus")

rank_abundance_plot(df_gemini1.0_pro_t1.0)+ggtitle("Gemini 1.0 Pro")

rank_abundance_plot(df_gemini1.5_flash_t1.0)+ggtitle("Gemini 1.5 Flash")

ICD converted responses

rank_abundance_plot(df_gpt3.5_icd)+ggtitle("ChatGPT 3.5 ICD")

rank_abundance_plot(df_gpt4.0_icd)+ggtitle("ChatGPT 4.0 ICD")

rank_abundance_plot(df_claude3_haiku_t1.0_icd)+ggtitle("Claude3 Haiku ICD")

rank_abundance_plot(df_claude3_opus_t1.0_icd)+ggtitle("Claude3 Opus ICD")

rank_abundance_plot(df_gemini1.0_pro_t1.0_icd)+ggtitle("Gemini 1.0 Pro ICD")

Combined model data

Original responses

multi_ranked_abundance_plot(df_gpt3.5, df_gpt4.0, df_claude3_haiku_t1.0, 
                            df_claude3_opus_t1.0, df_gemini1.0_pro_t1.0)+
  ggtitle("Combined model rank abundance")

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
Please use `linewidth` instead.

ICD converted response

multi_ranked_abundance_plot(df_gpt3.5_icd, df_gpt4.0_icd, df_claude3_haiku_t1.0_icd, 
                            df_claude3_opus_t1.0_icd, df_gemini1.0_pro_t1.0_icd)+
  ggtitle("Combined model ICD rank abundance")

Top diagnoses plots

Individual model data

Original responses

n_diag <- 25
top_diagnosis_plot(df_gpt3.5, n_diag = n_diag)+ggtitle("ChatGPT 3.5")

top_diagnosis_plot(df_gpt4.0, n_diag = n_diag)+ggtitle("ChatGPT 4.0")

top_diagnosis_plot(df_claude3_haiku_t0.1, n_diag = n_diag)+ggtitle("Claude3 Haiku t0.1")

top_diagnosis_plot(df_claude3_haiku_t1.0, n_diag = n_diag)+ggtitle("Claude3 Haiku t1.0")

top_diagnosis_plot(df_claude3_opus_t1.0, n_diag = n_diag)+ggtitle("Claude3 Opus t1.0")

top_diagnosis_plot(df_gemini1.0_pro_t1.0, n_diag = n_diag)+ggtitle("Gemini 1.0 Pro")

top_diagnosis_plot(df_gemini1.5_flash_t1.0, n_diag = n_diag)+ggtitle("Gemini 1.5 Flash")

ICD converted responses

custom_labeler <- function(x, wrap_width=33) {
    x %>%
        str_replace("___.+$", "") %>%
        str_wrap(width = wrap_width)
}

custom_text_formatting <- list(
  theme(axis.text = element_text(size = 7, lineheight = 0.7), 
          strip.text = element_text(size = 7),
          axis.title = element_text(size = 9)),
  tidytext::scale_x_reordered(labels = ~custom_labeler(., wrap_width = 45))
)

n_diag <- 25
top_diagnosis_plot(df_gpt3.5_icd, n_diag = n_diag) + custom_text_formatting + ggtitle("ChatGPT 3.5 ICD") 

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.

top_diagnosis_plot(df_gpt4.0_icd, n_diag = n_diag)+ custom_text_formatting+ggtitle("ChatGPT 4.0 ICD")

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.

top_diagnosis_plot(df_claude3_haiku_t1.0_icd, n_diag = n_diag)+ custom_text_formatting+ggtitle("Claude3 Haiku t1.0 ICD")

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.

top_diagnosis_plot(df_claude3_opus_t1.0_icd, n_diag = n_diag)+ custom_text_formatting+ggtitle("Claude3 Opus t1.0 ICD")

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.

top_diagnosis_plot(df_gemini1.0_pro_t1.0_icd, n_diag = n_diag)+ custom_text_formatting+ggtitle("Gemini 1.0 Pro ICD")

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.

Combined model data

Original responses

multi_top_diagnosis_plot(distribution_vis = "points", wrap_width=45, n_diag = 25,
                         df_gpt3.5, df_gpt4.0, df_claude3_haiku_t1.0, 
                         df_claude3_opus_t1.0, df_gemini1.0_pro_t1.0)

Warning: Duplicated aesthetics after name standardisation: size

ICD converted responses

multi_top_diagnosis_plot(errorbar_type = "range", wrap_width = 33, 
                         df_gpt3.5_icd, df_gpt4.0_icd, df_claude3_haiku_t1.0_icd, 
                         df_claude3_opus_t1.0_icd, df_gemini1.0_pro_t1.0_icd)

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.

Cumulative top frequency plots

Individual model data

Original responses

cumulative_frequency_plot(df_gpt3.5)$plot+ggtitle("GPT3")

cumulative_frequency_plot(df_gpt4.0)$plot+ggtitle("GPT4")

cumulative_frequency_plot(df_claude3_haiku_t1.0)$plot+ggtitle("Claude3 Haiku")

cumulative_frequency_plot(df_claude3_opus_t1.0)$plot+ggtitle("Claude3 Haiku")

cumulative_frequency_plot(df_gemini1.0_pro_t1.0)$plot+ggtitle("Gemini Pro 1.0")

# Example of cumulative frequency values
cumulative_frequency_plot(df_gpt4.0)$data

ICD converted responses

cumulative_frequency_plot(df_gpt3.5_icd)$plot+ggtitle("GPT3 ICD")

cumulative_frequency_plot(df_gpt4.0_icd)$plot+ggtitle("GPT4 ICD")

cumulative_frequency_plot(df_claude3_haiku_t1.0_icd)$plot+ggtitle("Claude3 Haiku ICD")

cumulative_frequency_plot(df_claude3_opus_t1.0_icd)$plot+ggtitle("Claude3 Haiku ICD")

cumulative_frequency_plot(df_gemini1.0_pro_t1.0_icd)$plot+ggtitle("Gemini Pro 1.0 ICD")

Combined model data

Original responses

multi_cumulative_frequency_plot(
    n_diagnoses = 25, distribution_vis = "std_error",
  df_gpt3.5, df_gpt4.0, df_claude3_haiku_t1.0, 
                                df_claude3_opus_t1.0, df_gemini1.0_pro_t1.0)

ICD converted responses

multi_cumulative_frequency_plot(
  n_diagnoses = 25, distribution_vis = "std_error",
  df_gpt3.5_icd, df_gpt4.0_icd, df_claude3_haiku_t1.0_icd, 
                                df_claude3_opus_t1.0_icd, df_gemini1.0_pro_t1.0_icd) +
  ggtitle("ICD converted responses")

Diagnosis rank table

Individual model data

Original responses

diagnosis_rank_table(df_gpt3.5, "mast |mastoc|anaphylaxis") %>% mutate(diagnosis = substr(diagnosis, 1, 60))

diagnosis_rank_table(df_gpt4.0, "mast |mastoc|anaphylaxis") %>% mutate(diagnosis = substr(diagnosis, 1, 60)) 

# diagnosis_rank_table(df_claude3_haiku_t0_1, "mast |mastoc|anaphylaxis") %>% mutate(diagnosis = substr(diagnosis, 1, 60)) 
diagnosis_rank_table(df_claude3_haiku_t1.0, "mast |mastoc|anaphylaxis") %>% mutate(diagnosis = substr(diagnosis, 1, 60)) 

diagnosis_rank_table(df_claude3_opus_t1.0, "mast |mastoc|anaphylaxis") %>% mutate(diagnosis = substr(diagnosis, 1, 60)) 

diagnosis_rank_table(df_gemini1.0_pro_t1.0, "mast |mastoc|anaphylaxis") %>% mutate(diagnosis = substr(diagnosis, 1, 60)) 

diagnosis_rank_table(df_gemini1.5_flash_t1.0, "mast |mastoc|anaphylaxis") %>% mutate(diagnosis = substr(diagnosis, 1, 60)) 

ICD converted responses

diagnosis_rank_table(df_gpt3.5_icd, "mast |mastoc|anaphylaxis") %>% mutate(diagnosis = substr(diagnosis, 1, 60)) 

diagnosis_rank_table(df_gpt4.0_icd, "mast |mastoc|anaphylaxis") %>% mutate(diagnosis = substr(diagnosis, 1, 60)) 

diagnosis_rank_table(df_claude3_haiku_t1.0_icd, "mast |mastoc|anaphylaxis") %>% mutate(diagnosis = substr(diagnosis, 1, 60)) 

diagnosis_rank_table(df_claude3_opus_t1.0_icd, "mast |mastoc|anaphylaxis") %>% mutate(diagnosis = substr(diagnosis, 1, 60)) 

diagnosis_rank_table(df_gemini1.0_pro_t1.0_icd, "mast |mastoc|anaphylaxis") %>% mutate(diagnosis = substr(diagnosis, 1, 60)) 

multi_diagnosis_rank_table <- function(search_pattern, ...){
  listN(...) %>% 
  lapply(., diagnosis_rank_table, pattern = search_pattern) %>%
  mapply(function(x,y) {mutate(x, model=y)}, ., names(.), SIMPLIFY = F) %>% 
  bind_rows() %>% 
  pivot_longer(contains(c("mcas","kawasaki","sle","migraine")), names_to = "criteria", values_to = "rank") %>% 
  filter(grepl("mcas", criteria)) %>% 
  format_models() %>% 
  format_criteria() %>% 
  pivot_wider(names_from = "model", values_from = "rank",names_prefix = "model_") %>% 
  rowwise() %>%
  mutate(mean_rank = round(mean(c_across(contains("model_")), na.rm=T)), 0) %>%
  mutate(ranks = paste(c_across(contains("model_")), collapse = ", ")) %>%
  mutate(output = str_glue("{mean_rank}\n[{ranks}]")) %>%
  select(Diagnosis = diagnosis, criteria, output) %>%
  pivot_wider(names_from = "criteria", values_from = "output")
}

rank_table <- multi_diagnosis_rank_table(search_pattern = "T78\\.2 |D47\\.02 |D89\\.41 |D89\\.49 |D89\\.4 ",
                                         df_gpt3.5_icd, df_gpt4.0_icd, df_claude3_haiku_t1.0_icd, df_claude3_opus_t1.0_icd, df_gemini1.0_pro_t1.0_icd)
rank_table  

rank_table %>% 
  flextable() %>% 
  width(width = 30) %>% 
  align(j = 2:3, align = "center", part = "all")

Diagnosis	MCAS - Consortium	MCAS - Alternative
T78.2 Anaphylactic shock, unspecified	1 [1, 1, 1, 1, 1]	151 [216, 87, 96, 186, 168]
D47.02 Systemic mastocytosis	8 [22, 8, 6, 2, 2]	56 [92, 76, 24, 45, 41]
D89.41 Monoclonal mast cell activation syndrome	40 [128, 23, 28, 11, 11]	61 [141, 68, 23, 39, 36]
D89.49 Other mast cell activation disorder	152 [307, 63, 96, 137, 155]	661 [1156, 557, 176, 383, 1031]
D89.4 Mast cell activation syndrome and related disorders	452 [725, 178, NA, NA, NA]	1380 [NA, 870, 1845, 1424, NA]

NA

Diversity

Individual model data

calculate_diversity_ci <- function(df, b = 100, seed = 1234){
  require(entropart)
  set.seed(seed)
  
  df %>% 
  count(criteria, diagnosis, sort = T) %>% 
  group_by(criteria) %>% 
  nest() %>% 
  mutate(boot = map(data, function(x){
    x <- deframe(x)
    x <- entropart::EntropyCI(entropart::Shannon, 
                              Simulations=b, 
                              Ns = x, q=1, 
                              Correction = "None")
    x
  })) %>% 
  mutate(shannon = map_dbl(data, 
                           ~log(entropart::Diversity(deframe(.), q=1, Correction="None")))) %>% # Calculate shannon with entropart
  mutate(boot = map(boot, ~quantile(., c(0.025,0.5,0.975)))) %>% 
  unnest_wider(boot)
}

plot_diversity_ci <- function(df){
  plt <- df %>% 
    format_criteria() %>% 
    ggplot(aes(x = criteria, y = shannon))+
    geom_point(size = 1)+
    geom_errorbar(aes(ymin = `2.5%`, ymax = `97.5%`), width = 0.3)+
    theme_bw() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    labs(x = "", y = "Shannon diversity") 
  return(plt)
}

Original responses

calculate_diversity_ci(df_gpt3.5) %>% plot_diversity_ci() + ggtitle("GPT3")

Loading required package: entropart

=============================================================================================================================================================================================================
=============================================================================================================================================================================================================
=============================================================================================================================================================================================================
=============================================================================================================================================================================================================
=============================================================================================================================================================================================================
=============================================================================================================================================================================================================

calculate_diversity_ci(df_gpt4.0) %>% plot_diversity_ci() + ggtitle("GPT4")

=============================================================================================================================================================================================================
=============================================================================================================================================================================================================
=============================================================================================================================================================================================================
=============================================================================================================================================================================================================
=============================================================================================================================================================================================================
=============================================================================================================================================================================================================

calculate_diversity_ci(df_claude3_haiku_t1.0) %>% plot_diversity_ci() + ggtitle("Claude3 Haiku")

=============================================================================================================================================================================================================
=============================================================================================================================================================================================================
=============================================================================================================================================================================================================
=============================================================================================================================================================================================================
=============================================================================================================================================================================================================
=============================================================================================================================================================================================================

calculate_diversity_ci(df_claude3_opus_t1.0) %>% plot_diversity_ci() + ggtitle("Claude3 Opus")

=============================================================================================================================================================================================================
=============================================================================================================================================================================================================
=============================================================================================================================================================================================================
=============================================================================================================================================================================================================
=============================================================================================================================================================================================================
=============================================================================================================================================================================================================

calculate_diversity_ci(df_gemini1.0_pro_t1.0) %>% plot_diversity_ci() + ggtitle("Gemini t1.0")

=============================================================================================================================================================================================================
=============================================================================================================================================================================================================
=============================================================================================================================================================================================================
=============================================================================================================================================================================================================
=============================================================================================================================================================================================================
=============================================================================================================================================================================================================

Combined model data

Original responses

calculate_shannon <- function(df){
  table(df$criteria, df$diagnosis) %>% 
    vegan::diversity()
}

tibble(model = c("ChatGPT 3.5", "ChatGPT 4.0", "Claude3 Haiku", "Claude3 Opus", "Gemini 1.0 Pro"),
       data = list(df_gpt3.5, df_gpt4.0, df_claude3_haiku_t1.0, df_claude3_opus_t1.0, df_gemini1.0_pro_t1.0)) %>% 
  mutate(shannon = map(data, calculate_shannon)) %>% 
  select(model, shannon) %>% 
  unnest_wider(shannon) %>% 
  pivot_longer(-model, names_to = "criteria", values_to = "shannon") %>% 
  format_criteria() %>% 
  ggplot(aes(x = criteria, y = shannon, color = model))+
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.3, aes(group = 1))+
  # geom_boxplot(aes(group = criteria))+
  # geom_jitter(size= 0.5)+
  geom_point(size = 1, position = position_dodge(width = 0.75))+
  theme_bw() +
  scale_color_brewer(palette = "Dark2")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  ggpubr::geom_pwc(aes(group = criteria), 
                   method = "wilcox.test",
                   label = "p.signif",
                   p.adjust.method = "BH",
                   hide.ns = T,
                   vjust = 0.5
                   )+
  labs(x = "", y = "Shannon Diversity", color = "")

                   # remove.bracket = T)
                   # symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 1), 
                   #  symbols = c("****", "***", "**", "*", "ns")))
  # ggpubr::stat_compare_means(comparisons = list(c("SLE - EULAR-ACR", "MCAS - Alternative")))

ICD converted responses

calculate_shannon <- function(df){
  table(df$criteria, df$diagnosis) %>% 
    vegan::diversity()
}

tibble(model = c("ChatGPT 3.5", "ChatGPT 4.0", "Claude3 Haiku", "Claude3 Opus", "Gemini 1.0 Pro"),
       data = list(df_gpt3.5_icd, df_gpt4.0_icd, df_claude3_haiku_t1.0_icd, df_claude3_opus_t1.0_icd, df_gemini1.0_pro_t1.0_icd)) %>% 
  mutate(shannon = map(data, calculate_shannon)) %>% 
  select(model, shannon) %>% 
  unnest_wider(shannon) %>% 
  pivot_longer(-model, names_to = "criteria", values_to = "shannon") %>% 
  format_criteria() %>% 
  ggplot(aes(x = criteria, y = shannon, color = model))+
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.3, aes(group = 1))+
  # geom_boxplot(aes(group = criteria))+
  # geom_jitter(size= 0.5)+
  geom_point(size = 1, position = position_dodge(width = 0.75))+
  theme_bw() +
  scale_color_brewer(palette = "Dark2")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  ggpubr::geom_pwc(aes(group = criteria), 
                   method = "wilcox.test",
                   label = "p.signif",
                   p.adjust.method = "BH",
                   hide.ns = T,
                   vjust = 0.5
                   )+
  labs(x = "", y = "Shannon Diversity", color = "")

                   # remove.bracket = T)
                   # symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 1), 
                   #  symbols = c("****", "***", "**", "*", "ns")))
  # ggpubr::stat_compare_means(comparisons = list(c("SLE - EULAR-ACR", "MCAS - Alternative")))

Permutation testing

# Run outside of notebook
# Includes all 10,000 GPT iterations and 1000 bootstrap permutations 
# Contained in script source(here("scripts/diversity_analysis/permute_null_diversity_difference.R"))
library(here)
source(here("utils/data_processing.R"))
library(tidyverse)
library(vegan)
library(broom)
library(here)
library(furrr)
library(future.apply)

model <- "gpt4"
p <- 1000
i <- 10000

print("### Reading data")
read_path <- sprintf("data/processed_diagnoses/diagnoses_%s.csv.gz", model)
df <- read_csv(here(read_path))

print("### Calculating permutation")
perm_out <- difference_permutation_test(df, metric = "precision", permutations = p, gpt_iterations = i)

print("### Writing data")
write_path <- sprintf("data/diversity_analysis/precision_permutation_test_%s_p%s_i%s.RDS", model, p, i)
saveRDS(perm_out, here(write_path))

```r
diagnosis_pca_plot(df_gpt3) + ggtitle(\GPT3\)

<!-- rnb-source-end -->

<!-- rnb-plot-begin -->

<img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZgAAAFQCAMAAAC1Xpi+AAADAFBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29wcHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGCgoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OUlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWmpqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnKysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7////isF19AAAACXBIWXMAAA7DAAAOwwHHb6hkAAAcHUlEQVR4nO2deWAUVbaHKxiCiewMyBMkbIMoMAo4mnkuSSCCgkDGIM4omzqIQRQcdYILREAQJQiyibIkPhkUiWyK7ATDlmB08MloIMCEyN5ZOvueM3VrvdVddXtJJV0J5/ujq+v+bnVV95eu6u5UncsBYkk4X28Aog+KsSgoxqKgGIuCYiwKirEoZoqxv/HH1r+7d1Exub+KI3R8+DjAWE5mOUD8gObBE34zcaWNFRPFHGvnPyQmJozrdQ2ImBfi4t6bFBTwb9gRFxcXyb3I3/4LlnERH85s1bPIvLU2VswTc7FV1xQy/brZw0DEHCYzh7mJQricEzLoEF4NsIn7p2lrbbSYJ2YSd0i881JgliIG2t8jTCQxl7n5/O1v3LumrbXRYpqYytaD6VlJTFHTIcKsJKYqu4y/3cAdNGutjRfTxJzhXqVnBTHVp0dxi4VZeVfGf0KYNjLg1Sqz1tp4MU1MGreITLqST18L5E9lHPfXaiFVxVwacHPAqxVmrbXxYpqYHO4FMvkkLm6h3wLpU1ncR6lSqorh+YybYdZaGy/mHfw7/K905z/iO+YwHUpiyu01ZHJ3P9PW2mgxT8wUbrd4Z5GxmM+4NDJ54C7T1tpoMU9MTtuux8l0d6CxmHTuH/xtZuAk09baaDHxm//BW24Y9sYbEU3e72EoBp7wG79q7i1tLpi31saKmb+V2V7+002/eywFPthiKKZsVo9m3caeN3GljRX8ddmioBiLgmIsCoqxKCjGoqAYi4JiLAqKsSgoxqKgGIuCYiyKWWJOfWVI4ibjjE+/SGTFm9jpRlb6lYv0S2bMTr9gPynvn/JRk8UkfLzXiF07DCOSbt7Fir9hp5tZ6d7tzHQbM97DXnjzt6x0JzPdtXm3Ybb7ZbPF/GAYVZSwFqywVbLiYnaazdyofGZqL2ClNeyFbaWstJyZVtiMT0apQjEoBsU4gGJQDIrRUg9ijk0dO/MawMmp41dXG0wEUAxN3Yu5+sSpmoTZUPX0T6X/SNKfiKAYmroXkzQL4NqT8MNrAIdi9SciGjH57z29sFDdDhSjxQwxpfzre3Am7FwKkDlFfwJQlJCQMPdoicK12zmOu8Mm3L9y6dJv5y8x+C3jAis+z0wzM/TbC8QNySthkcuMi9kL2wpYaSE7tRUZZ9PcE8OTPOFX+HI1v1cbqz/hb8LDwycnZSt8IJxOvoS/tz/0qed8wKRBbwsbYstmwU59tPDVF90UUzB3WibAt8v4N8fz+hMRelcWI4h5AyAnNMel97rhrXhy24h3ZRXT15LTwH+YAXA0Vn8iQotZL4j5AmDdevaj1x3lw8htIxbz3eti9wlnq2ft15+I0GIqH+K9PMyv/D3fXTnW2MWsHcHzFMCvL00mbx3diYDmU1n19rhvSLAAxejh+2/+kpiO3C/CtDuXViZeh+GEUbvCjN6k15sDbuo66RJ/L2Sq0Np5AcCj/NvzxpDlZPY/o9q2/8slIUExEkwxfrPJ5McmXFrlq/rFF4zaFYiYsru6Ld+3+vY7CrViwlNSdr7elP/cXt3/gaSDd4UKCYqRYIq5X7iC7M0HnN4WxdrZmmowgoiZ04V8yLO1XKEVM5rc+9LvJPybOwWwi8sj8yhGgilmYdN0fnLbEi6tindji2ozYAuXDf6pQ4fB2cgOLUJPAGn33xLs130j/3JGd2k54rTDQ/FialoI+yvY9LWOmJqer0HmEl7s1iByxTqKkWGKWfPIXICTN2aIYkIeTtn4OyIm5NMs6BO2Lyl8gCjmlk0nn2xWDIPCDh4fd7PDjy28mEzuR2U2ZEwaoYMiBh57jNzufLvTfGEWxUiwxay7E2B25BVBzHcBNnKdJi9mDv/0484AbGgtiuHtneXSU5rm8oeLTonah+LFHOAuK7Mh0tXqqpjnQ8jttD6tE4RZFCPBFpPT9DT0Wy+KWX47kBozvJhv+DtlOxeMaSuJ2cG/WFx6vPiSC1VNEvk7J4XH4MWckd4x2ad1dmXSOwZgu3ggQzESbDEwdN7pALsoZvEdfNNRIoZPC/v3mZP0pSRmnyBmc4cqgvA5oMJut4sfCHgx1UHLhLt/uVdPTK/XIG0judPuY3KLYiRciFnTf/5wEMXsbsYfPhZLYrYH8PutjRoxp/1+BrgU8r32ocinsphbSMGtrFZv6YjZzC+1tU05/wj+35J5FCPhQky2/63xkpiae4anbe7Jf6glYg5zn2QlBje7SomBqB47DoT3cqhqQsQU9+687MCK7l0LtGIGpaXtiw2IBshr89fjR4b+Xvj3D4qRcCEGHvLPkcTAlZGtwjb51whi4J2ObaPS+/ajxRRHd2496pzDQwnf/IteviOw53PkbeP4zf9eYS+XOrz9zWMyhQTFSLDEaMlZVQSwPtik9RqBYiT0xSTpNJa2e/nyz31nm7ReI1CMhK6YjzfpdT0WEtQ9psyk9RpQPZTcohgDMRcfYm5gHbLiA3KLYoxOX9oV+voCHxAbNU349oNiDM8rq/gpLS3lSBqDlH2prPgwO92v324T145i8IQ/LSgGxbgCxdCgGBTjChRDg2JQjCtQDI2FxKxLrTCipNAwIqmtlBUXMtOCbFZaYWemecy4PI+5sK2IlRYz0xJbmWFWOt1sMUeNryzI9/aaBJ58Zmq3sdKSXHbKvgyDvbAFLsNwD9yV0VhoV4ZiaFAMinEFiqFBMSjGFSiGBsWgGFdcJ2IOP/e4yMjHHmcQxU5HjtZpfUa4ahLFeCXms7Fni3WazaAmc9JKQDEEz8VUhdXlcIQ1ESUohuC5mKznmA9YW177GcUQPBfznynMB6wtr59AMQTvxawZ2KJNBDn9t680KGukeE1PmOMCTh0dLr6m76IYGa/FvB+0JO27yTccocQ8mEL4xaG/c0eHi6+pu4BiZLwW0yOW3Eb8jRITqfsYzh0dLr6m7gKKkfFaTLsnye2vx12Kce6ovfiavg4bUIyM12JmcncuOiHUalHEhAongDqOIuncUXvxNX0dNqAYGe8P/t9OuJVr92K+08HfaQhpp47ai6/p67ABxcjU6uPy6ZiWYQa7MvVSa+eO2ouv6euwAcXIeCvm5GjhJ8r1nE1fjHKptU5H7cXX9HXYgGJkvBWT5beVTFY2r3Rx8Nfp6HDxNXUXUIyM17uySUFzdn/3QctYfg81lhzzT8oHf8caRs4dHS6+pu4CipHxWkzlorB2re/+mP+41Vc45veWD/6cw2+czh0dLr6m76IYmYb+W1kyKYWAJeQJlhKTEVkMWEJe5Le/MR+wtrx60gMxc6JG8GI8LCGvpfGIqQ6rcG40jepBZR7tyqJ4MewS8oTrQgwkjv7/3Lril6fiwWMx7BLyOdHR0a8k243IyzWMeHJteaw4h53aWKk9m53meL5w8gsTRcZPmMhgAjOdOF6vMXqf8Ixe8kgMu4R8XgzPoQIj7HmGEUltdlacl89Kc22stCCHnebWYmEb+0kx03yb8ZPK80yMZyXkHWhEuzIVC5xXRsR4VkLecTtQjBZaTPWGqbH0fws8FONhCXmH7UAxWmgx48gvAqfUDL/5W0PML8JPNU+rGYqxhpivBTEPqBmKsYaYs4KYaDVDMdYQAzN4L90vqhmKsYgYSJ6/ht4OFGMVMQ6gGBSDYhxBMSgGxWhBMSjGFSiGBsU0CDFnRt967xZqHsVYQ0zRbeRnl/1qA4qxhpjdwu9hE9UGFGMNMZsEMY+pDSjGGmIuCGJWqg0oxhpiYFMnjptEDYqLYiwiBkp+ukrPohiriHEAxbghptToxMic88zzJjMus9IrzPRqhk0zX0NtFIoRxGQ8OopVv6qeGPSeulEohogpeTCL2aO+mKN+ZEMxRMzWpcwO9UbVEPUuiuHFfLyV2aH+GKbcM13M2iNFRhTYDSOS2gpYsZ2Z5tlYaVEuM83J/dAyYpSNyje7hHz88UojSosMI5LaylhxITvNZqWVdmaal79SEuNumSu1OUy8tqzrO+oy/Ff8JoFCGaWxpNONd25QF5QTtQqWts8wZaPKzB50oUHuylaJYtwuc6U264pZ1Nzvn2Q69o98p71/8TsKDglVBUvbpw53ZQ1ZjNtlrtRmXTF3T/jTCDIdG0FuK1rEgkNCVcHS9kExNIoYt8tcscVkcDveD8gD+UWHjspnPimhq2Bp+6AYGkWM22Wu1GY9MXNal5/m1oH0ohcsbn8OtAldBUvbB8XQKGLcLnOlNuuJuX0iwB3kG4lwYOea7FOWkxK6Cpa2D4qhUcWAm2WumLuyE9zKM2f+5n9NOrDvj2pbJC0qJ3QVLLUPAcXQyGLcLnNFNQ8ZJ0z+531VTIz4dvpIPn7kcD9Ii8oJXQVL7UNAMTSyGLfLXFHN028jpRh+4XYpYmqCBVe3hSoH9qafg0NCVcHS9kExNMquzO0yV2rzuVaDv0peFfxglbLMEfE8mHlNLsovepvF4lJqQlXB0vZBMTSKGLfLXFHNv47pEnj7m+Q/OtIyU28VdndZTZbIL/oD/cSlqEStgqXtg2JoNAd/34JiaFCMZcV8Hs/sUH/g/2No7AXZg5ibVm8kvqHcdRaTs37RcTK9uNyrx26QYmBvaOwC3zNxXJmyUU5iTtwc1J4jlX8OevceaphioOzYXgP2bDNKBBJ3sNKd7DRxNz27L5PaKCcxEcOLarY2++h6E2OMVc4ra5nC33zS6iqKkbCKmO7byMbcH1mDYkSsIia6yzf8kzzbYkoiihGwipjSR7m7+UlKJz8UI2AVMQBZ/xIW2hHH3CAjUAwNnlR+HYgpXr0bIOO+uV4OWYtiaEwUc7lf80386iJv6sl+rjTXXW1/inoTM7m3OFbGuS4vMDeIfozrrrY/Rb2JCV4t3Vnai7lBFNdfbX+KehMTJBcI2B/E3CCK67C2v0q9iblLvnZmWR/mBlEotf0vDxw48NkDNsQErrzoIGZGj1xheu3W6XoS9FBq+xd/9dVX7x4rNaKowDAiqa2YFeezUxsrLc1jprnMuIS9sI35pAq9fspFjpdhlPYLXvHj5R+XdurJfIPTYG1/Q8z8HlMS25zjuMBXipjbo3kMrO1vhLnf/GvOHThTrdvbAKztb4S5YjyS4giKoTFRTMEzrQMGn2VuCxMUQ2OimBeaz4wL/oP3bxoUQ2OimM78t8UU7gxzY1igGBoTxfgdBqjkUpgbwwLF0JgoRnCCYlRQDIqhcRbzYGRkpHCjf82OK1AMjYliRqswN8gIFEOD//NHMa5AMTQoBsW4AsXQoBgU44qGIabk68/V+jAoxjJiMntyHPelPIdiLCPmceGSfPm6OhRjGTHdBDE/SXMoxjJi7hfE5EhzKMYyYrYTL5PkORRjGTFwcGzkR8oDohjriNGAYlCMloYtpjjhrRiev78Ww+DvL7HSmOns9GV2/C7r34bXrZjLYevSfMyh52cZb3aDEbMm2W5EXq5hxJNry9NrHv8vkzasNkw8arjZNvaTYqZ5+k9ZIOclad1Wfcc8XMvtMYVtKw2jBvOOMVvMML3G+mb3YsMIxfgSFOOMKKYjJw5I0Z1Lo8f4UMf9UEcYcUDuwB0Gm6srfvgOZZxjxVgCinFGEuM3m0x+bMK/btQYH+qIIOoIIw7IHXgxYQtcbDvfofLV33QCFOOMJOZ+ofztmw/wYqgxPtQRQdQRRhyQOziLIRUmajRnzxuaQzHOSGIWNk3nJ7ct4dLoMT7UEUHUEUa0KB24w3dzXATkR3dpOeI0gH/q0GHgvyXYr/tGgLORHVqEngDSoYpLG0nOdFzVtlzuCihGr1kSs+aRuQAnb8zg0ugxPtQRQdQRRrQoHbjDVaHzqmFQ2MHj427OBv+QT7PA/5ZNJ59sVgx9wvYlhQ8A0oEX889A/r00aIrSFVCMXrMsZt2dALMjr3Bp9Bgf6ogg6ggjWpQO4q4spWkuQHWnRPCfw7f587LPcuk1cWcANrQWdmW8mIIbt8DVG1KVroBi9JplMTlNT0O/9bwYeowPdUQQgjjCCNDjiagdRDHxYun3d8H/Gz713wFg59KhbOeCMW0VMfDYeFh5m9oVUIxesywGhs47HWDnxdBjfCgjglAjjJCHssvjiahDhohiNneoIlSDP/lk7b9PEFPYv8+cpC9VMV+0rQybr3YFFKPXrIhZ03/+cODFUGN8qON+UCOM0FBDhohiTvv9DHAp5HuNmO0B/E5royqmKGi9/3m1K6AYvWZFTLb/rfGCGHWMD2pEEHWEERqqAy9m0ORciOqx40B4rwqNmMPcJ1mJwc2ukg5EDDzefjCfyl0Bxeg1K2LgIf8cQYw6xgc17oc6wggN1YEXk9BuFBRHd2496hxoxMA7HdtGpfftRzoIYhK5T/lU7gooRq8ZfytjPYIHNMqf/Xd+aBhdt2L+7H4FyLrjnd2G0XUr5uijl2q5RbWmZttw44IU160Y+P7JCEL44AgG4aGsNCKMnYaz44diGRV1r18xEnj6Ui1BMTQoBsW4AsXQoBgU4woUQ1MvYpLJ73lSGX/diQCKoakPMRmRxUoZf92JCIqhqQcxc6JGFCtl/HUnIiiGpl52ZVHFShl/3Qm/GampqYtTK4woKTSMSGorZcWFzLQgm5VW2JlpHjMuZy9sK2KlJezUVmaYlcqF4t0TI5Xx151gbX8zcartr8Ou6OgfJTFSGX/dCf/2u3Dhwqrvq4woKzKMSGorZ8WF7DSblVbZmWlePiutZC9sK2alpcyU9ZTLPXvHSGX8dScieIyhqbdjjFTGX3cigmJo6k2MXMZfdyKAYmjwmz+KcQWKoUExKMYVKIYGxaAYVzQwMVWnL6MYC4o53ovjRlxEMawFfSGmqge5SGMyimEt6AsxvwhXz/RCMawFfSHmvCCmP4phLeiTY0wYEfMeimEt6BMx157pfPsHeSiGtSB+j9GAYlBMXYuxJ+01YpthQtjKjPcYpZnCmlGMKzGfRcxbUK+Me468rCjGhZhTf67VENHe8NESQDEuxSzcw9ysuqCKXLiLYlyIiXUqvFT3oBhAMbXGx2J6zyC3k1p6oCmds5OJTtFFuWxjyFRhtvMCuhFQjIAnYqY3dy7BaIwkRqfooly2USNGbgQUI+CBmDeCvmP1dUQSo4NctlEjRm4EFCPgvph5gcLJh1I9xXunA4zjLsAVbrfcBNvvCgyOU3oQMWmtllU57cqUso20GKURGpiY+LQaI8qLDSOS2ipYcTEzLcqeJYtZItZ+k+spzroToJv/57AxsFRuOh/wj+NxXLI8y4v5qe374CxGKdsYMkao8d9hAd3IM4xfs62EtV1lzLTcVmmYVZotZu2hAiPseYYRSW12VpyXz0pzba9LYnq3HtKZ/EHL9RSP+uVcDHgiGqYMU5r2cacAtp2TZ9O5ox1iQEeMUrYxRCxHxy2gG3mG8Gu2sZ8U8znl24yfVF7d1PbPj9f5DWPeXNYvHPNmzmfFc9nprIipwp32zabEthxA7r3zzCN/CAo8W9Vm88b7EvpA36Ugl1gseejGyOV2ZTada9c+ShGjV7aR2pVpij02qF2ZIObX0M+YPxvWG9tGbnj8xRdnZPql+2WAUmIR/v3WgBaJ8mw6tzCZ2yWL0SvbSIlRG6EhihlxwWW/eqIqYlm//jug67ge/DFfKrGY9DYfjBkhz5KD/7jflznuyqiyjaoYqhEaoJhqS1R/E1n4uZ9/Hky8gX9t5RKLSdyik4nBs+VZIuZKq3ccxVBlG1UxVCM0QDEVo0x6MBNYtqPvH/it4kg5X7nE4ofdmwXPqJRnhe8xSwPPOIihyjaqYqhGQDG1YtmOunvshirG67FF6L6uMRxcRKAuxTwCDVSMl2OLaPq6xnBwEYE6FFNARmdokGK8HFtE01emxuAflMUuhoWpOzFVk8h4HA1SjHdji9B9s8d17PhX/tucPKCI9CuX0pw6dJg0uEgR+VqYyV0B/xWdg0KzXunYfqHwKMvuf7xuiAr/P/L4DVKMd2OLUH1r7rkn6WBI/xp5QBHpVy61OeTTLGlwEVVMt+SDnQNj0qP9hIOU8I7B05ccxHg1tgjVN+mG8/xnU//98oAi0q9cajMZbESsyK+KSQCY1qMGLnHCgK8oRkIjxquxRai+K3uSO79fLg8oIv3KpTaTbyeOYvYBvD2Y/EcSxdBoxHg1tgjVd4VgoPcSpQq/+CuX2kw+alNiMmQxESjGEa0Yb8YWofoe8M8CuNh0ryxG+pVLbabFHOI/L6AYQ7RivBlbhOpbc/efDh+5jxzlJTHir1xqMxEjDS7SftjJpG4oxhCtGG/GFqH6gu2pjjcLn4ulXZn0K5fSTMRIg4vs6t3khlUoxhAf/lZWpvcSohgJIqZyhOt+9cXSnYBiCMI/ygYZr7C+eTYDUAxBEJMQzRg/oj6pTphAJihGPhljy6NDnMf7GDyINRrI4FD2QCReDVMy/F1h8DcUg5f6aUExKMYVvhCT+dnnhv8oE2jMYo5NHTvzmjVr+2/gOK7rOdbCjVjM1SdO1STMtmRt/3LhN9Eo1sKNWEzSLIBrT1qytv/PgpiurIUbsZjSQoCDM13U9i9YunRp7JEiIwrshhFJbQWsOM8wvSiIuY+1cA4rLMrJZcfM1JbPSvNdPOVC4yWnuSeGJ3nCry5q+18bOXLk1IO5RuTkGEYktbFj4/RpImYja+FsVpib7SJmpl5vtbAwI3K3tn/B3GmZrmr7E3xw8K9aM3r8t6xlG/OurGK6UIsca/s74msx370udsfa/g74WszaETxPYW1/J3wtxm1QDA2KQTGuQDE0KAbFuCJhfaoRR5INI5LuOcqKk9npXlaamsRMDzDjFPbCew6x0iPsdM8xw+yo2WJ+TjAkfp1xlpCw9NnlrHhtPCudPZmVJqxhpn+PYcZrWWH8s/NZ8TrmU/7g2Y+MQ/n7h1livObEwHTvF/4kvBZrnhTj/bLVAzd5v/Cxgeddd0IxXoFi2KCYuqP8Qrn3C+dfrMWar+XUYuELRd4vW3aB+VFTxOdiEH1QjEXxhRh3T/AwIJmczOftwgqe9K3tar15wj4Q4/YJHvpkRJIxhr1cWMGTvrVdrVdP2Adi3D7BQ5c5USPIqNzeLaziSd/artarJ+wDMW6d4MGADP7s9cIynvSt7Wq9esK+Ofi7PsGDAXmFvF5YxpO+JqzW8ydcz2I8OsFDd2HxFfJ8YQc86StSi9V684R98I5x+wQPA8gr5PXCMp70re1qvXrCPhDj9gkeBpBXyOuFZTzpW9vVevWEfSDG7RM8DCCvkNcLK3jSt5ar9eoJ4zd/i4JiLAqKsSgoxqKgGIuCYiwKirEojV3MaI7j/LpMJef2lb7eP6jb81fE9hg3iwj7jEYv5sGUlENLmj8LUNKv+4p9q3v3E6qqHObY53D6nkYvhlSvhrdaA8QGk5MvrjVfDbAnqhmK8TGimLiAqprmUlnn3fz7JW4aivExREzV8eBBcI4uAQ1pKMbHjBYuOu93FjQloFGMzyEH/5SMarWss41UmEMxPkc8xvBUB4plnSPvI7coxtcoYuCVTqQA97mb5pIZFONrVDGFPbssP7A8uKfwPQbF+BpVDBRMuyOw5/NC3XoUg3gLirEoKMaioBiLgmIsCoqxKCjGoqAYi4JiLAqKsSgoxqL8F7G/V2WbqJisAAAAAElFTkSuQmCC" />

<!-- rnb-plot-end -->

<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuZGlhZ25vc2lzX3BjYV9wbG90KGRmX2dwdDQpICsgZ2d0aXRsZShcXEdQVDRcXClcbmBgYFxuYGBgIn0= -->

```r
```r
diagnosis_pca_plot(df_gpt4) + ggtitle(\GPT4\)

<!-- rnb-source-end -->

<!-- rnb-plot-begin -->

<img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZgAAAFQCAMAAAC1Xpi+AAADAFBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29wcHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGCgoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OUlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWmpqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnKysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7////isF19AAAACXBIWXMAAA7DAAAOwwHHb6hkAAAcLUlEQVR4nO2deWAURbrAKzgEEjkiCPLkCAIqAllEWY3rkYCBcCyQJ4jKgscKIgiCT1wQUCQcRgjCBlAuOfahKxKBRRAQJICBDRgR14ghBIQoZw5iCCSQ49uumulzeibVPZ1Mk3y/P6a6q6uqj99Md3VNdxUBxJYQf28Aog+KsSkoxqagGJuCYmwKirEpVorJn/zHkFsfnHeFTi8hlGa9DgEMJSKLWLJLz39v4UqrKxaK+XdjR8+JEyPJXReBinklPv69EcGBP8HW+Pj4GDJW+HQKeYZssW6l1RbrxJxp2DqFhl/U6QVUTDKdSSbPs4WLSIqY7v8JiuHAOjEjyDfOiVeDsiQx0OQBFshiTjXAXwwPlokpCXlcOesSU1i7J5uVxJQ90jcJxXBgmZhMMkE5y8SUZQwg89msJGbmref2oBgOLBOTSubRoDWtfcWJtTJCniljS0Uxh2pvABTDg2VicskrNFgWHz83IM5VK4v/8KBrqUtM4Z1CXQDF8GDdxb/pn1wTvzh/McnKhS4xcY6/r149iYxfXWDZaqsr1okZTXY4J+Z5FjNVvNXMtGy11RXrxOQ2an2IhjuCPIth4KmMBwvv/PfcflOfyZOjas1pi2J8x8q2suzXHrr51idS4P2NKMZnsHXZpqAYm4JibAqKsSkoxqagGJuCYmwKirEpKMamoBibgmJsilVijn3ujfWJXhd7Yt06U9kS15vKtv5TU9k+N7mRn+ofkwMWi1m9dKcXtm73ttQjmzebyrZ9i6lsWzeYyrbT5EZu0D0mO16zWsx33pZevW6q0N9/N5Xt+hVT2a7mmMoG5jayJLtEL7oUxbiBYrhBMb6DYgyCYrhBMdygGN9BMQbxi5jf33th7mXlUhTjhj/EXO5ACOlUyKbz8yhnL+SZ4fRpU9kunNHd54qo/mKWsoctVwhTqVF/eckPjOg+x8QuVH8xE5mYyQC5EbkWFW+UqauM56n+YtYyMZ8CrFxrUemGudbHeJ7qL6akh+ClVynAe3ssKt04KMaFqlZWtjl+S7kQxqEYj/j1PsYlphk5ysI2JLWYpOpm9xQvMak9TTXlvptbjzgrTIWPYbEt4gD+LPw864azLgR+GdCoydNn2RIU48KrmIDpNDhci6SWTPhVN7uneAkqpvjeOxbtWn5Ph8tqMd1SUra9WXu08EPt8mjSnnsj2BJrxPw6qM/AJyuJbu+V62xE1Yp5JIwGUx51+1lo7gTLyzyWT8XEtqKVvOwGi9ViBtGpzwLS4CdyDGA7uUTnLRFzvVu68WK4mZGgE1m1YubWpjt49wKSWiq4yR54y30bSQ44Dkb3gRMxTetHHAEa79gYGtBmnXBqGNWqQb8MTVGCmPL6zi5P1n+hI6a83RtwaoEgdlNwMZ23RMzOd42Xwk9ZD53IqhWzovcMgLS6x51iwnulrLuViglfkwUdI3cldbvPKeb29WlD6lyB7pF7Dg27TXOUBDGnyGFpNnxwKqWpJAaeeIJ+bnun+Ww2a4mYNf80XooB9LaxisWs7AwwPeY8E7M3MJu+pymIiRW+6PGZAJ+EOMUI9k6Q9JTaecK3qXmiuihBzG5yTpoNd73LKYt5OZx+jusYsprNohgX3sXk1s6AsLVOMYvuAdrHjCCGvlhWvC1ucCOXmK0A+SR9lfOQs7NIojCRxsoQxGS6fjE5GTqnMtcvBmCz80JmqZgV99e/JYruSydXJxMxzm2M1GZwS6ipTConbSEGomdlBOY7xczvIEQdoGKEpZe7dIxN+swlZhcTs6FpKYXVA67n5+eXiXtYFryQTT79oJ6Yu96A1HV0ovFS+mmlmDnBC1L3jrxpv0LMYymUo5r07gk1lUnFpIdtrGoxK7rM7gtOMTvqCAdgvkvM5kDhvLVOJSYj4EeAs+HfqouiezjxdtrhVlbDqTpiNgi5Nt1yTSjB8SWdt1JM22n0M2q4QkyMbgnuCTWVScWkh220XMyqQyXuzJbE5DharnKJKX+gb+qGdkKllopJJsuyEkPrXFCIgYFtt+7udpfmHxy6h1fat1i4e3Gb1gVqMd1TU3dNCxwFcOmWZw7tj77zKl3SW2dzKqAwRxuz0imm8RD6+fOhCsW4J1RXJpX1SoE+OhtRPN5VllViPtpf6E6sJAZ6OHJdYuB8/4aR6x3lTAzMbNZoYHqnMKWYK6NahAw4qSmfnawLX+sQ1O4l+rPR3vk/yM5yB/s2uW3wKbYkWmdzKiA/RxuzzCnmLdJ53hF2NyiJiWCVwt80G+meUF2ZVNYrBfrobMTv41wLq+JUpiZ3SSHA2lCL1usJSy/+Xz7XkjQe+7vbxX+SNoNbQnVlUlmv9LCNVXSNSdKJLGr82rkfO023aL2esLq6nDGxQaSHU5lcdXRPqK5MKuuVHraxasQsXa+X9N/hwW0mFlu0Xg+URRvP40lM2qAiGqwl2fpipKqjTkJ1ZVJZrwR/ijnTo8ii4o2y+H3jeTyJyQrYRIMP6pVUcPHXSaipTComwZ9iYHvEm3F+YNrAcZ6bQz3i8VQ2Ijh2x973G0wTzlBD6TU/Tbz4a9tk3RNqKpOKSfCrGLj+A9uB/SmpZti711S2lORsM7vgUUzJvMjGIV2XCtWtTuya3168+JNSdXr3hJrKpHJSEHPNfSPwgT93qr6tbKp7HIpxp+rFBLr/3YNi3HEXs+FDUwXx0jPU/R8ZFOOOu5hL3a6aKomPr15dRxK1kSjGHZ3//HdFvru0klg4+n8vQ+8WlzVrRDHu6D0lc3Vfxa+5/svcy7HJJQAZdSZqVohi3PHH40uTHf9RR6MYd/wh5uodnSa/r1wvinHHLw/89RbuQ1v9IkejGHf8ISaDtRCMkKNRjDv+EPMFE/OoHI1i3PGHmGP4i6kYv1xjRuM1pkL887T/lknzsFbmnRr0GoYEinEDxXCDYrhBMb6DYgyCYrhBMdygmAp5vV+/fksA0sY8u7xMGzBQjEEsEjO04Pr1Uih94YeivyVpAicoxiDWiCkayoLv3gD4ZpomcIJiDGKNmF+Gjh0yMw+2JQCcGq0JnKAYg1gjJn12bmnCu/DZcoALQzUBwMX+/fuP2eOtB7HcXFMdj+XkmMpmcm252aay5ZnbyDz9tWWPNSSGcv4p+HKh8Bt5WRMAFCQkJEzTe3FJfieowPh7RAJ5eaayFVwylc39xSU+ck3lKsjWPSYGX1w6fgwgeyh8NwngwDRN4ARPZQax5lR2+IULZUsWQelzJ8re/loTOEExBrGouvz58OfmC7v786sjPyrXBgwUw4c0LDve+XNTBWL+fhu56yvnJIrhpvLFbGX/8GexaRTDTeWLeZmJcXZ2g2K4qXwxo5mYf7BpFFO2dQ7fm5sz3jb3xud0jjTzDtJN2Um9tHJ28VPjxVzt+66pp/GtZdsr7ICvaEk673VuV40XM9Pt3SD/MGE3C6QDUePF9DHxYnll8K3mTVgUY6oI60l7Qz2PYkwVYT0oRl8Mb++JcnTkcBbReqacB+C3WkHsrcmhNFHdzp/IGcUlcueK6jQoRlcMd++JcrSumHn1Aj6m4dA/Col2Ph1wADRLFJ0rqtOgGF0x3L0nytG6Yro+91A/Gg6NYiuqPw00SxSdK6rToBhdMdy9J3oXc5xsnRNIe0J3HnRoJnU27lqi7FxRnQbF6Irh7j1RjtYTExtyLYOsBNdBL5jfROoh0rVE2bmiOg2K0b/48/aeKEfribnneYAOPcF1YSe1dkn5XEuUnSuq06AYj9Vlrt4TvZ7KjpAPMjOHOy66LuxfD2xU6MoqLlF2riinoaAYPTHcvScqonsOY8H/zJHFOAdfIx+K149c8p0rq7hE2bminIZtAorREcPde6IievzddKOPku2SmPJQ5uruCOnCXtvVb5a8RNG5ojoNitE9lXH3nihHn2z4+Of7loQ+Virl2U/Yoyezap0RD/ot85255CWKzhXVaVCMrhju3hMV0T8PbhV0zxT69IQrz5iW7HSXVWuBeNAfDXPmUiyRO1dUp0Ex2FamwrZiepsqwnp+/Jt6vsaLeeGEqTIsZ+0y9XyNF5PW86ypQizmuwhNX341Xgx8/2R0FBfdI/jSaYnkSBM98pxmI1EMNzfks8sVgmIMUlViVuzL90LeJW9LPZKbaypbnsls2aay5Ztb26Vs3WOS+6rFYvAXYxA8lXGDYrhBMb6DYgyCYrhBMdygGN9BMQZBMdygGG5QjO+gGIOgGG5QDDcoxndQjEFQDDcohhsU4zsoxiAohhsUw41VYsqP8nSSsOVzc50rcAy3+PV/yrUbiWIAzkSPM9cZiXVM6K59vA3FAAzQvgHrB05HaSJQDOQ9ZaoUi3lZ85NBMXBqlKlSLGbKYfU8ikExXrCDmGbEealpQ1KV/VfIfVrIvWdoEBOQZMgucF+sQkhQTLRvQ1FMi8ldO+8QDc8sqmDV+twAYgKm0+BwLeG4KfqvkHu7kHvP0CAmEMRExlWwbiFByYRfdRaYFXPktuAmhD4FuMfcb+gGEPMIe7VryqOCGEX/FXJvF3LvGRrEBO5i6IrKVf1ueTRnVkxU38LyTXU+rM5i5tZOF4K7F5BUZf8Vcm8Xcu8ZaqQEJLkrIVHw+6hWDfplADgORvcBx8bQgDbrAE7ENK0fcQRoglKS2p++Yruk0TUxKZgX0yBF+FjW8EI1FrOi9wyAtLrHSaqy/wq5twu59ww1UgKSXBoxqwy6R+45NOy2HHCEr8kCx+3r04bUuQIdI3cldbsPaAJBzMdBwiZ0Hy0lBfNi2vxL+Ch/JKbcgJgbZMQlSczKzgDTY86TVGX/FXJvF3LvGWqkBM5TWUrtPICy5ongiBXiHILsEyS9PD4T4JMQdioTxBTU3QgXbjooJQXzYka12iKIPVF/dCK3mBtlxCVJTG7tDAhbK4hR9l8h93ZBcfaeAcq+MuQETjGrnK81vwuOLcJSx1aAfJIOxdviBjeSxMATz8IHd8tJwbyYoj+TrkKQ0jyAW8yNMuKSJAaiZ2UE5gtilP1XSL1dKHrPAGVfGXJ3GE4xG5qWUsrAQWvWjl1MzOUuHWOTPpPFfNqoJHK2nBR8uY/J+p5t0NZ43v2XRlyy+fgxRyUxK7rM7guCGEX/FXKfForeM5QousNwiskI+BHgbPi3KjGbA4WT1jpZTGHwWsdpOSkVozlA1owfo8uNMuLSD5KYHEfLVUyM3H+ForcLufcMJYoEgpjuI/NgYNutu7vddV0lJpksy0oMrXOBJqBi4MkmjwtLxaRUjPYA8Y64dGX5DqFm+PAM/j+ZpBGXKDfEqQx6OHKZGLn/CkWfFnLvGUoUCQQxqxsPgCujWoQMOAkqMTCzWaOB6Z3CaAImJpGsEZaKScH8qexcWL31wi1XzM3tuP/Hu1FGXLqx28pGtnfWS062eoV3VTfKiEunR5oqxWImmxQTutw1kXAX97pukBGXCgeYKsVihp1Rz/OKCRa/+l8Hm1qvjcXAC9tMFWMpyYM0Ebxi7v3ANbGwo6kV21nMlZG9n/QzfZ69pNlIXjGT2uax8GLL8aYOhp3FCEeBp5Z9NtNU5TzvNEeaYreN5L7zDwtdfPjc4YTm7Sr6K0gfe4vhwqbPlV2dVo8QEvR6oanVVB8xp39w/3Z7p9If+Cs/uTvT9GAr1UTM1f6EhO40lq/Sxfg0Ak41EfMmba9saeysUbliCv4aEvi4D53hVRMxD7KWZPe//b1RuWJeqfdWfOgfzP9oqomYHkxMmqF8lSumRQJACsk0tQ5KNRHzMfXykLHvZ+WKCUgW0pIUU+ugVBMxsLprm+HnjeWrXDHMCYoxA4rhpXqJeSwmJoZ96A9GUBEoxiC8YgbJmFoPijEIPlTODYrhBsX4DooxCIrhBsVwg2J8B8UYBMVwg2K4QTG+g2IMgmK4QTHcoBjfQTEGqSoxH9noxaXfL1WcRof8HFPZCk2+XVV5Ly6pWHWoxAuFRd6WeiQ/31S2osumshXmmMpWYm4ji7OLdaPFB2HxVCZRPU9lKMYgKIYbFMMNivEdFGMQFMMNiuEGxfgOijEIiuEGxXCDYnwHxRgExXCDYrhBMb6DYgyCYrhBMdygGN9BMQZBMdygGG5QjO+gGIOgGG5QDDcoxndQjEFQDDcohhsU4zsoxiAohpsbUszr/fr1WyIN5aMOGCjGIBaJGVpw/XqpOJSPOnCCYgxijZiioSxwDeWjDpygGINYI+aXoWOHzMwTh/JRB8KBOHr06EJ8DcMY1ryGkT47tzThXXEoH3UAcO7+++9/cXc2YgHn3Qb28cj2UaPYCBvnnxKH8lEH+IsxlcuaX8zxYwDZQ8WhfNSBE7zGGMSaa8zhFy6ULVkkDuWjDpygGINYVF3+fPhz869IQ/moAwaKMQje+XODYrhBMb6DYgyCYrhBMdygGN9BMQZBMdygGG5QjO+gGIOgGG5QDDcoxndQjEFQDDcohhsU4zscYr7sF2UPxusfypoqZvMwk8M4W87uXrrjKdZUMb2vWrQq35n+tV5sjRVj0ZosYMNyvdiaKqaPRWuygI1L9WJrtphm5CgL25BUgOIp993cesRZFvFbraDLbGLF/fVvidrjVoQybcVkF0AxXYM+KEZEISZgOg0O1xIOW/G9dyzatfyeDkzIvHoBH9NwTvCC1L0jb9IOEK5KWzGRcVAy4VdPS1GMiELMI2E0mPKoICa2Va4wmd1gMY3p+txD/WjYdhr9jBquKUGVVqTcw2jVV6gYL6AYEYWYubXTheDuBSS1vP4iFrX+C+HjONk6J/CSMNF4CI37+ZC6AGXanGHNmj1zEcCxMTSgzTqhHn5vUGi8IvpgdJ+uhESVktRCOgD8KXIeHItbBEdkvd6syVxWCooRUYhZ0XsGQFrd4yT1FDksp4gNuZZBVgoTb5HO846UawtQpC1/4IGkPeFdysFx+/q0IXWunA7826F4sk+ODl+TVRoxq0wl5o59e1oETUwfFcAuUihGRClmZWfhRiLmPEndTc7JKe55HqBDTzr15XMtSeOxmr1VpE266TRAluNrcAiGT5D0XeQYwL9OytGxwK4xKjGrAca1LYezzmG+UYyIUkxu7QwIWyuIyXT9CnIyAI6QDzIzhzsuOpNnTGwQySYSCaHHFkCR9oN2dOLOReDYCpBP0q/2qBuzKF8RvQXcxewCeOdxoQqBYtQoxUD0rIzAfEFMWfBCFvX0gwATCeNDSBtUROPWkmwaXM/Pz2cXeEXaxcxA+wXsaAtiAH6ael/9RDmaVrUVYo6LYqJqpJiKx4+RxKzoMrsvCGJg4u30F5LVcCqUhw6jC++OgKyATXTqg3qazZbT7nZkAZypvVMUk/SOsGBwPzlaKeYbob6gI+ZDacuKPp8tEvv2bFO8YyTxnP2WvoZRMRWPuCSJyXG0XMXEXGnfYuHuxW1aF8B+wlqvZtU6AyOCY3fsfb/BNE35ctryrg8l73+YXuVdYsi8tMTQ6XI0FdN9ZJ4gBpr0SUu6Q0dMgrhh56Pe31mlbB89rmpHXDJwKoMejlwmBgpf6xDU7iXhpzCmJTthZdVaACXzIhuHdF3qVi+T0kL2X5rdxurFrlPZ39vUCZ1UIkdTMasbD6BitrevddMSb6eyKdss2n9+Xv2GBXa7xlQ1xXpXAFlML7cvQKWTPIMFNV2MLrIYP2zUd1NZYB8x0RatyQI+XiNOoRiAFz029lY15TGnxEmnmPaT6OeIBoc8ZNAhneTTQKcJW2wmDx/DZlvEKSPBhmIudJ9ftRUgT2x8Yom0YQox4+tpW7S94RKj04QtNpOrxIiRYEMxcG3T3DhjTJ9uLP1LjpaThWDWDK+plmXKGyaLmRy818j+usToIDaTq8SIkWBHMcYx/PjSgWa3pxh6fEkSMyuI3UmdiGlaP+IIPCjc5g0jvwmV+h1ilNiI7ZqlYlIbLix1O5VJzeRKMVIk1FAx8Nsf66wyI2YBeZdNdYzcldTtPni7M8Adjn/CuqAiMcrViC3OCmJ+aDQH3MVIzeThg1MpTeOUkVBTxUDRs+Slq4bFtA/p2YJ+ocvjhdPcJyFwICD3TOBTo2B0HynK1YgtzqaTA00ngo4YqZk83Nn8R+KUkVBjxQAsqNWT98kASUz9f19q/iydKt4WN7hRCJTesmHdw6s7QqcEKcrViC3OppPGTQZKYqRWcEUzueJUpmo7r7Fi4MuQtj+d3bBbd7+1uMS8IeQimwAud+kYm/RZCMCTY8dOOhWQHnBcjnI2Youz6WTuPrJdFCO1gsvN5EoxciTUZDHw0111hcPwh4scSRXV5b/elg2bA/MA1gkWVoR12Qqth7UFKcrViC3O0ov/sDuLtacyuZlcIUYRCTVazPVj7As6giOpQkx+y0GQTJZlJYbWuQC/BjguwfM3CcdWjHI1YouzVMz5hjO1YhTN5LIYRSTUbDG7mZg/cCRV3vnvIP+Emc0aDUzvFAbQSci9mtA/Q8UoVyO2a5bdxyQEZWrEKJrJZTGKSKjZYo4yMb04kmJbmalCTb8f8wQVs50jqR8eqP72LRbUTDFF70U/5f6UrQ5/yTK1Bl9Y/g8W1Ewx3Em/78V/z2MNB7o5X0dBMd45PDhafNGse4S5F9QijSTuO9b1LA6K4aZ6Pr6EYgyCYrhxifl1bSLnux4uUAwvPonZIFS0W6cbyYdiePFFTFlzegtk6PYGxfDii5gTrNGgqZF8KIYXX8RcZWK6GsnnLzH76EHVHdGnWg7s839UTKKRfH4SczxG2E3dEX2q58A+ZWuefD7JUD7/iIkd2O+KhxF9cGAfJ/46lQ0UdlN3RB9pYB8KijGINWJ0R/SRBva52L9//zF78ryQm+ttqUdyckxlM7m23GxT2fLMbWSe/tqyOQb2cY3oQ8XojugjDexTkJCQMK3CF5dMkJdnKlvBJVPZ8nNMZSvMNZXL9xeXqBjdEX1wYB8n/jyV6Y7ogwP7OPGnGP0RfXBgH4bN7/zXHvTCvv3elnokKclUtv37TGXbt9NUtoPmNvLAVwd0o60W8+Nqb6xc5XWxJyZMMJVt1Uemss0cYSrb6hWmci16caFuvHhtsEpMpTB+bMVprGPNn6pybT/d/6PX5ShGAsVwg2JsSo7JepI5Cs5U5dqu/XbN63Jbi6nJoBibYj8xfH/MWU/llk4xtGe2E8P5x5zlVG7pFGN7ZjcxvH/MWU7llg6G98xuYnj/mLOcyi2dYWjP7CnG+x9zlULlls4wtGc2EsP9x1ylULmlMwztmY3EuOD7Y85yKrd0hqE9s6eYiv+Ys5zKLZ1haM/sKabiP+asp3JLpxjaM/uJQRgoxqagGJuCYmwKirEpKMamoBibUt3FDCKEBLQaUyBMFr3ZJfiOl8874ycaezO56qn2Yh5LSflmQb0XAa6GtVm8a3n7MPaYZjKp0scJTFDtxcTQz6khANNC6biAF+stB/hqYB0U42ecYuIDS8vrucYF3CH8XuLHoRg/Q8WUHgrtDieVYwhCKorxM4PY6+FhJ0A1hiCK8Tv04p9yvEweFzD7OP1EMf7GeY0RKAtyjgsY8zD9RDH+RhIDrzenIzievJmNRIVi/I0s5nK7Vot2Lwptx+5jUIy/kcVAwbgOQe1eZgOfohjELCjGpqAYm4JibAqKsSkoxqagGJuCYmwKirEpKMamoBib8l8WvF9qZtV33gAAAABJRU5ErkJggg==" />

<!-- rnb-plot-end -->

<!-- rnb-chunk-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxucGVybXV0YXRpb25fdGVzdF9wbG90KGRpdmVyc2l0eV9wZXJtdXRhdGlvbl9yZXN1bHRzKVxuYGBgIn0= -->

```r
permutation_test_plot(diversity_permutation_results)

Similarity

Jaccard/Bray curtis

diagnosis_similarity_heatmap(df_gpt3.5, method = "bray")

diagnosis_similarity_heatmap(df_gpt4.0, method = "bray")

diagnosis_similarity_heatmap(df_claude3_haiku_t1.0, method = "bray")

diagnosis_similarity_heatmap(df_claude3_opus_t1.0, method = "bray")

diagnosis_similarity_heatmap(df_gemini1.0_pro_t1.0, method = "bray")

diagnosis_similarity_heatmap(df_gpt3.5_icd, method = "bray")

diagnosis_similarity_heatmap(df_gpt4.0_icd, method = "bray")

diagnosis_similarity_heatmap(df_claude3_haiku_t1.0_icd, method = "bray")

diagnosis_similarity_heatmap(df_claude3_opus_t1.0_icd, method = "bray")

diagnosis_similarity_heatmap(df_gemini1.0_pro_t1.0_icd, method = "bray")

• Bray-Curtis similarity measures the similarity of a given diagnostic criteria’s set of alternative diagnoses along with their frequencies.
• This demonstrates that SLE criteria results in a very similar set and frequency of diagnoses, while the diagnoses associated with two MCAS criteria are as different from each other as they are from those generated by the criteria of other conditions.

PCA

diagnosis_pca_plot(df_gpt3.5) + ggtitle("GPT3")

diagnosis_pca_plot(df_gpt4.0) + ggtitle("GPT4")

diagnosis_pca_plot(df_claude3_haiku_t1.0) + ggtitle("Claude Haiku")

diagnosis_pca_plot(df_claude3_opus_t1.0) + ggtitle("Claude Opus")

diagnosis_pca_plot(df_gemini1.0_pro_t1.0) + ggtitle("Gemini")

diagnosis_pca_plot <- function(df){
  df <- format_criteria(df)
  diagnosis_table <- table(df$criteria, df$diagnosis)
  diagnosis_pca <- as.data.frame(diagnosis_table) %>% 
    pivot_wider(names_from = "Var2", values_from = "Freq", values_fill = 0) %>% 
    column_to_rownames("Var1") %>% 
    prcomp()
  
  as.data.frame(diagnosis_pca$x) %>% 
    rownames_to_column("criteria") %>% 
    ggplot(aes(x = PC1, y = PC2, label = criteria))+
    geom_point()+
    ggrepel::geom_label_repel() +
    theme_bw()
}

multi_edge_density_plot <- function(...){
  df_list <- listN(...)
  df_list <- lapply(df_list, ungroup)
  
  # Process data
  df_combined <- df_list %>% 
    mapply(function(x,y) {mutate(x, model=y)}, ., names(.), SIMPLIFY = F) %>% 
    bind_rows() %>% 
    nest(.by = c("model", "criteria")) %>% 
    mutate(edge_density = map_dbl(data, ~edge_density(graph_from_data_frame(.)))) %>% 
    select(-data)

df <- listN(df_gpt3.5_icd, df_gpt4.0_icd, df_claude3_haiku_t1.0_icd, df_claude3_opus_t1.0_icd, df_gemini1.0_pro_t1.0_icd) %>% 
  mapply(function(x,y) {mutate(x, model=y)}, ., names(.), SIMPLIFY = F) %>% 
  bind_rows() %>% 
  count(model, criteria, diagnosis) %>% 
  pivot_wider(names_from = "diagnosis", values_from = "n", values_fill = 0) %>% 
  unite(id, model, criteria, sep = "__") %>% 
  column_to_rownames("id") %>% 
  prcomp(scale. = F)

as.data.frame(df$x) %>% 
    rownames_to_column("id") %>% 
  separate(id, into = c("model", "criteria"), sep = "__") %>% 
  format_criteria() %>% 
  format_models() %>% 
  ggplot(aes(x = PC1, y = PC2, color = criteria))+
    geom_point()+
    # ggrepel::geom_label_repel() +
    theme_bw() +
  scale_color_brewer(palette = "Dark2")

Precision

• Precision represents how similar each iteration of a 10-point differential diagnosis is with all other differential diagnoses from the same set of criteria.
• I.e. how reproducible the 10-point differential diagnosis is for each criteria
• Measured by obtaining the Bray-Curtis similarity values between all iterations within a criteria

# Script for calculating all Bray-Curtis similarity values within a criteria
# Found in source(here("scripts/diversity_analysis/calculate_precision.R"))
library(here)
source(here("utils/data_processing.R"))

models <- c("gpt3", "gpt4", "gpt4_icd")

for (m in models){
  print(sprintf("READING IN DATA FOR: %s", m))
  read_path <- sprintf("data/processed_diagnoses/diagnoses_%s.csv.gz", m)
  df <- read_csv(here(read_path))
  
  print(sprintf("CALCULATING PRECISION FOR: %s", m))
  df <- calculate_precision(df)
  
  print(sprintf("WRITING PRECISION DATA FOR: %s", m))
  out_path <- sprintf("data/diversity_analysis/diagnosis_precision_%s.csv.gz", m)
  write_csv(df, here(out_path))
}

summarise_precision <- function(path){
  df <- vroom::vroom(here(path))
  df %>% 
    mutate(distance = 1-distance) %>% # Convert to similarity
    summarise(mean_cl_normal(distance), .by = criteria)
}



df_precision <- data.frame(model = c("gpt3.5", "gpt4.0", "claude3_haiku_t1.0", "claude3_opus_t1.0", "gemini1.0_pro_t1.0")) %>% 
  mutate(path = map_chr(model, ~str_glue("data/diversity_analysis/diagnosis_precision_{.}_icd.csv.gz"))) %>% mutate(data = map(path, summarise_precision)) %>% 
  select(-path) %>% 
  unnest(data)

df_precision %>% 
  format_criteria() %>% 
  format_models() %>%
  ggplot(aes(x = criteria, y = y))+
  geom_point(aes(color = model), position = position_dodge(width = 0.75))+
  theme_bw()+
  theme(axis.text.x = element_text(angle= 45, hjust = 1))+
  labs(x="", y = "Average Bray-Curtis Similarity") +
  ggpubr::geom_pwc(method = "wilcox.test", p.adjust.method = "BH", hide.ns = T, label = "p.adj.signif", bracket.nudge.y = 0.3, vjust = 0.6, step.increase = 0.14, tip.length = 0.02) +
  ylim(c(0,1)) +
  labs(color=NULL)+
  scale_color_brewer(palette = "Dark2")

df_precision <- vroom::vroom(here("data/diversity_analysis/diagnosis_precision_gpt4.csv.gz"))

df_precision_summary <- df_precision %>% 
  mutate(distance = 1-distance) %>% # Convert to similarity
  summarise(mean_cl_normal(distance), .by = criteria)
df_precision_summary

calculate_boxplot_stats <- function(df){
  df %>% 
    nest(data = distance, .by = criteria) %>% 
    mutate(box = map(data, function(df){
      x <- boxplot.stats(df$distance)$stats
      x <- sort(x)
      names(x) <- c("min","lower","median","upper","max")
      return(x)
    })) %>% 
    select(-data) %>% 
    unnest_wider(box)
}



manual_box_plot <- function(df, width = 0.9){
  df %>% 
    format_criteria() %>% 
    ggplot(aes(
      x = criteria,
      ymin = min,
      lower = lower,
      middle = median,
      upper = upper,
      ymax = max
    ))+
    geom_boxplot(stat = "identity", width = width) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))+
    labs(x = "", y = "Bray-Curtis similarity")
}

df_precision_stats <- df_precision %>% 
  mutate(distance = 1-distance) %>% # Convert to similarity
  calculate_boxplot_stats()
df_precision_stats

# Run outside of notebook
# Tried multiple combinations of bootstrap permutations and gpt_iterations
# Because a single bray-curtis matrix of all 10,000 comparisons takes ~10 minutes
# Contained in following script:
# source(here(\scripts/diversity_analysis/permute_null_precision_difference.R\))
model <- \gpt4\
p <- 1000
i <- 10000

print(\### Reading data\)
read_path <- sprintf(\data/processed_diagnoses/diagnoses_%s.csv.gz\, model)
df <- read_csv(here(read_path))

print(\### Calculating permutation\)
perm_out <- difference_permutation_test(df, metric = \precision\, permutations = p, gpt_iterations = i)

print(\### Writing data\)
write_path <- sprintf(\data/diversity_analysis/precision_permutation_test_%s_p%s_i%s.RDS\, model, p, i)
saveRDS(perm_out, here())

Permutation testing

# Run outside of notebook
# Tried multiple combinations of bootstrap permutations and gpt_iterations
# Because a single bray-curtis matrix of all 10,000 comparisons takes ~10 minutes
# Contained in following script:
# source(here("scripts/diversity_analysis/permute_null_precision_difference.R"))
model <- "gpt4"
p <- 1000
i <- 10000

print("### Reading data")
read_path <- sprintf("data/processed_diagnoses/diagnoses_%s.csv.gz", model)
df <- read_csv(here(read_path))

print("### Calculating permutation")
perm_out <- difference_permutation_test(df, metric = "precision", permutations = p, gpt_iterations = i)

print("### Writing data")
write_path <- sprintf("data/diversity_analysis/precision_permutation_test_%s_p%s_i%s.RDS", model, p, i)
saveRDS(perm_out, here())

readRDS(here("data/diversity_analysis/precision_permutation_test_gpt4_p1000_i10000.RDS"))

readRDS(here("data/diversity_analysis/precision_permutation_test_gpt4_p1000_i10000.RDS")) %>% permutation_test_plot()

iNEXT

```r
ggsave(plot=plt_top,filename=here(\figures/3_Top_25_diagnoses.pdf\), width = 7.4, height = 6.5)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->



# Final figures

### 2_Top_diagnoses

<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuY3VzdG9tX2xhYmVsZXIgPC0gZnVuY3Rpb24oeCwgd3JhcF93aWR0aD0zMykge1xuICAgIHggJT4lXG4gICAgICAgIHN0cl9yZXBsYWNlKFwiX19fLiskXCIsIFwiXCIpICU+JVxuICAgICAgICBzdHJfd3JhcCh3aWR0aCA9IHdyYXBfd2lkdGgpXG59XG5cbnBsdF90b3AgPC0gdG9wX2RpYWdub3Npc19wbG90KGRmX2dwdDQsIG5fZGlhZyA9IDI1KSArIFxuICB0aGVtZShheGlzLnRleHQgPSBlbGVtZW50X3RleHQoc2l6ZSA9IDUsIGxpbmVoZWlnaHQgPSAwLjcpLCBcbiAgICAgICAgc3RyaXAudGV4dCA9IGVsZW1lbnRfdGV4dChzaXplID0gNyksXG4gICAgICAgIGF4aXMudGl0bGUgPSBlbGVtZW50X3RleHQoc2l6ZSA9IDkpKSArIFxuICB0aWR5dGV4dDo6c2NhbGVfeF9yZW9yZGVyZWQobGFiZWxzID0gY3VzdG9tX2xhYmVsZXIpXG5wbHRfdG9wXG5gYGAifQ== -->

```r
custom_labeler <- function(x, wrap_width=33) {
    x %>%
        str_replace("___.+$", "") %>%
        str_wrap(width = wrap_width)
}

plt_top <- top_diagnosis_plot(df_gpt4, n_diag = 25) + 
  theme(axis.text = element_text(size = 5, lineheight = 0.7), 
        strip.text = element_text(size = 7),
        axis.title = element_text(size = 9)) + 
  tidytext::scale_x_reordered(labels = custom_labeler)
plt_top

ggsave(plot=plt_top,filename=here("figures/3_Top_25_diagnoses.pdf"), width = 7.4, height = 6.5)

3A_Rank_abundance

```r
plt_div <- df_diversity_ci %>% 
  format_criteria() %>% 
  ggplot(aes(x = criteria, y = shannon))+
  geom_point(size = 1)+
  geom_errorbar(aes(ymin = `2.5%`, ymax = `97.5%`), width = 0.3)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = \\, y = \Shannon diversity\) 
plt_div

<!-- rnb-source-end -->

<!-- rnb-plot-begin -->

<img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAASAAAAEgCAMAAAAjXV6yAAAC61BMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycpKSkqKiorKyssLCwuLi4vLy8wMDAxMTEzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJERERFRUVGRkZHR0dISEhKSkpLS0tNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29wcHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGCgoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OUlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWmpqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnKysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7///8iNV5cAAAACXBIWXMAAA7DAAAOwwHHb6hkAAAWlklEQVR4nO2daUBUV5bHy1Fn7Gwd0ll60klPnDgznfQ4tpqZTpuxydg2aTdEkSBEWWIURBYZ48R048K4hGg0sV0iGjXu27ijVkQEF1BBRBFtiY2gVlFUFZSs1v04b6n3Csr77jm1GJF3/x/g8epw6/HjnXvPOfe+WwbCxZThUV9AZxcHBIgDAsQBAeKAALEA3Z4ZPrNOOiqdHLny/g9zQZ1NDEDO+DOtX30uHrVFFzemGV2nv5tOU3IK9XRHpSYjjNKS0xBWyakIoxQ/LmotDOh8OiEtFvGoKJWQE7Ncp7OLaMZ19Yj/hqMWYdRiakVYme8hjOxWhFGjyUk525bkOmAA2j/3f+Oy7OLRgSWEVE4SDhqys7Mz8u9RVGulnfWQzYwwqjc1IKzMNoRRnQVhZDM5aFeRCAPaMvysY/EC6WglIXcihIM7gwcPjjOaKTLRTvpmZcJYBe796G93ZwoMaG+GQCRMupe+EO6geNdp7mKKzswWAIWLR0VCb5Q/y3WaA1LUElnu/DKLkPJ7beOv3//kqOt0FwNUlzp8HeU0BhC5ODkyU+ikh5WRywlxqxXOXQiQMSgoqKfBYHgyKKiPx0soQHR1IUAlcXGxAh9Dn7i4aR4vcUCiDhtUveDxEgck6vaWLVvmvGwY+M2WLbs9XuKAVPk+immIAwLEAQHigABxQIA4IEAcECAOCBAHBIgDArT6ZANFtXW0sx6ymhFGdpMdYWWyIowstQgjq6mectaGKLlqaM2ZVoosNtpZD9XXIowaTU0IK3MDwshWhzBymFooZ5um+gyIuxggDgiQfgH9Ys4NRHs6BjT7zW6//tIMNqhfQISUzny9Z8gm4K31DIgQW0p3w9MTmbeRjgHdXf37v3sx/sjuQT9vYzSoX0DvdH95ynFxFZDNcJXRoH4BJZ10WbadYi2W0i+g30lf60YCDeoUUP3ChYaFoib9BGhQp4Duvv224W1Rg1YBDeoUkKBfeRolh4SELPc40i2gC7fIBZdUowh7S0ubx5FuARlSiDKPr/5mxINHRLeAmlpJk0uKzY2IKeGzLR2Omk+fPp11uoUii4121kP1ZoTRPVMjwspcjzCy1iGMGkzNlLONtILZyUqya/gc9f93ZW5t25J5HY5q+vXrN+GYSQe6TVnEuciw1/bMyFc+bn+n3Q7rcNRWVVW1/GwbRRY77ayHGmoRRk2mZoSVuQFhZKtDGDlMrZSzzZQ76NUFZNm/kc2vKT9XlBNiiuh4RHTbBwn60Vnyh5mktJfy87noO/eXLxUXcSpHkvQLqG96Ra9LJOsf1RPbJ47PckiLOF1HkvQLaE/Pbu+STw2LgQb1C4h8f+QeOXoMalC3gBzvVSGa0zEgErwR0ZyeAZ36j3kH8wQBDeoXkGcupiH9AkJKx4Bqd39GSlgTGpL0C6j8xZcMZOAblUCD+gX0+3FtPUjVO8OBBvUL6OmzpAchR4KABvUL6Ge5IqD9rwAN6hdQ7OC6HqSsz0dAg/oFZB3c0/BStxEOinl7fX2qiSKLlXbWQ3YzwshhciCszHaEkdWCMKo3NdKugrKIs4WcWrODNSuvd0DPJ1Cdx1P6dbFl7/zNG/NvgQ3qFxAhNUsH9RyyAWhQz4AIaV3/PE9WRVEB3dsd9dwToduABvULaOSTvUZsggZ5PQP6wwY7oj29AqKt7qBLp4Aoqzs0pFNAlNUdGtIpoCNuAQ3qFNBTgroZDN0NT7wJNKhTQIJ2v7Lb0fpdnx1Ag/oF9EsJzf/9Uj3hXrrZfidO/QJ6Jlf8mvesekJdutlhJ079Anr3XSFQtP/XEPU31aWbfCdOSRU/DRoyJOjvK5Sf3Ys41Z049b0EjzhWTk1e5f7PuBdxqjtx6nsRJ023O+7EqfNlwB5yL93kO3FS5V7E2ZV34vRndYd7EWcX3YlTkheAGq9JAhrUL6CtvXi5QxUNUO8x5dJe00CDjw+gixFjzsBWaEBt3S4hLuxxATSnd+9/6CG4w6u9eyezLdGAnK9uR1zY4wIoe/z4d8X+YtD48UvZlngX2/HyvEPfCQIu7fEAlOze6DeSbYkH1MMl4NIeD0B/LSwsTH8uKEH4doNtyVe5AvIG0PHlS4006w7SLyBT/+6vvd6jf9cZ5gMNKHxgJSE333ofaFC/gF6Uxq8TLwEN6h7Qi0CD+gUUPvCm4GL/Pg5oUL8bTZr69+jzeo9fmQBAOt6q1GlctuwYH+Yl8XqQKl4PAqTjelCAAXWtepCkwALqWvUgSQF2sS5VD5IUYEBdqh4kiT/1DIjHQYB4HASIx0GAeBwE6OHHQY175e98J05RlDhoSaz8/QfaibO5sLBw34Cnf31Y+A5skdE54qCCRBnQD7UTZ6V7ts9gYZt2ijjIklAqA3Iv5yQPFVBdZmbmCwKcV4TvAIBHB8jUXzlyzjr/vQxIXc75kD9O/Xy7O+iaXy15YYX/OPUL/Z8X9MQvlJ/3riAuQKKk5ZwN2dnZGfm0j7WvtUIffC/IZma/3lBdXX0t4Z20G8J3B9vUbEO8X50Fc1Em2lvVU0qug95Z9+zG9T+/ovycNWbM6PfGiHcp34lT0lNGMuo42Rbazky6gwK2E6d9XiL0OCPpzIBeOErmZpCb7fe0lwAFZCfOvcHBwT8Wupa3goOnsC07L6CQwX/Z92bdyp8CDfoGaPfQoYPFvvefhg59bAFV/MvCxn/9W8McoEHfAC13j09vsS07LyBhYCf1e05CDfoGqConJ2eqwTBwf04OsLCyEwNqrpQENOjHKGYuRxh1XkBbfsTrQapogF6bcOO2KKBB/QJ65jKiOT0DGgo97yxLp4Dy8vK2vjwvh+8G7NIDgHq4BTSoU0B46RuQHVw/pVtABR+eJ/VDDEGzIUQ6BZTb/bdXSOLT2UueWgk02CkAlUxf3gxbBRJQcJyQij2fRsjCAUCDq0/YKTJbaGc9VGdGGFlNVrbBsL59/1kI+H/ct2/fnWxLSy3mokw22tkED0DP5hJyziDkSsbnAECP/A5KCg0VARmGhYYCIUkg76AnhT97gbhT8qFO/yGQbwUFSVljUFAQsKVhIAH1XULIm+OFgz8NBhp85IBWZWZ++jODYVRmZmYZ2zKQgJY/+ekYwwHSsueJz4AGHzkgUfcLoZxaVCABOef1/slMQkZ1i4E+FqFTAHpkgWLpTbBBfQNCiAMCxAEB4oAAcUCAOCBAHBAgDqj9L7sWcfKdODXkWsTJd+LUkLKIk+/ESZe6iFPdiVMUB6TIvYhT3YnTHBERkXi8jiJzLe2sh2pNCCOLyYKwQhnVmjFG1Isye1YUKXIv4lR34rTNmTPn47x6imottLMeqjMjjGwmG8LKZEUYWWoRRlaTnXaW9ty8h9yLOPlOnFpyLeKEd+K8Nn02sDReVBcFBO3E+Ulw8G96GQxP/TY4eB37WroeIJoeADRlwIDXxCL6GwMGLGH/qk4BDXOvz/yU/as6BXQhJ+fQfxoMUTk5OdfZv6pTQJKuwh9tp29AOh3maeKAAHFAgDggQNmx0ymalkw766GUJIRR6tQ0hFVSCsIoeRrCKGUq7Wya74DoGvdHhNGy3yGMzvXDPLDwmzUIo7RY2Ibs6sf8RDUOiAMC9AMBWrwVYXQQehZN1F+mVyOs/seIMFr/Z4TRmektrJcDBajLigMCxAEB4oAABQZQ7WHMyIPSjVnp8D5Gpk3fQE+OBkoBAXQ1MiN0H2hVcvQ+aEOKIzZtHAuVt8sjv8rKhJ8qIUVramAjQAEBNLWEVE2EjJxR4bEgouJI4VZcVsA2ao5FbOAgqCx+QQRQwhN0OjZiPiOxCwigiGpSEz5nJZA/fvvZ6aS4Y0xEpWGiq6alZzE/t/xSuvTt/Cb2+91fe50UjYcIVcebGjNnab8eEEDfrCWff3l9ASMMLC4jpCH8DtkTEm9kIDLFrCJkY1LJF9Gsv+vCDOnb9ZGNrIsyjg7ZT2BCx5aR8oSaw5qv+w/IcbbA0ULaCGkd69A0OhImEFr/VXVs4YmY+DsaRkLXWxizamOy0MzeUC03E4yuhsqxdiTrb7d8dHNPZKVIiNGjHTlLKt6/mHCnZWSDlonfgEo/SE8dfUw8sr/PWIOeIxCyR8QUEdJ2XsNE6nrvxkRImLU2K5SMtieIf8+tMEYSVTMhk5A8kRArccmbRsi8UbeILVzz4Sp/ARVHlhByefymw4uOJx9gGR4eXUayFzIMXF2v5GWAkTMrpsBeNvkg6/2MI6tchBhyJp0l9/40aV3CXk0TPwFdHCsFQLfHFWxeUqxhUyL9ISs/DCuzhlVpN6V0vUtZhBSj+NQRkxi7jORuEwhNEMb4vJ0aFqVmcnBLM7kg3ELOc9uvaFgRPwEVl1VHbZSOjqZrW10LFzrLldMbBS9bPU/bTO16/xpDLep2NGpkhkHWyQohLR3Y2VyZEX3UObOQ1Q7xE5DQ9boIOYYxLlggJPAR+6FyRlGkWe16tfr6K/PbmjH9c4tKSKtQ3lRKbsabyeX/TtwKfH6dny6WoxCqYMaJ18bGSSMye8TdCXW9juSMNshIcJ3KScK4ZY0RCGl2QFVRy7bvixMGgqKpIUAY4GcfJHS9IiFHYq6GgZxfSF6mLdkI7noFQq2AkeA6ZJdIKD+E9TjiuvfyiETICSWRPgNyd73VUWuTN2pYKfkFk5CahBjZXa98D7GMXK4jEro7zcho6NLVqNMyIUg+A2rX9So9NU1KfnEtvELTpl0SAmWgIiGGkeI6u2L2TGds1dIkOHylRAieXfPdxdxd711G5qTmF6wMHU5CSo67DiRC2lJcx/jHI9pGxrGj80VCOTuuMppyyTdAUj0G1/VC+QXOqDRM2aTYkazV3QkqKEK4ztXYmvJhgpNWpsIlGh8BueoxQNcrC8gvsEYKISdhzEG0TrqLcJ0168iJjGhwMzuXfAGk1mNQhJj5BcrofmWDQsiRysgdqhIzCMJ1zn7b8JG9cFgC/HiuKF8AuesxzK7XJUZ+oZRXnSwjcj0pPvyGTMihOVxKyhkuujvTdUwilq1rSeUXpezLVuQLoHb1GMTaX6KZXyjlVUuqQ9uI1ERfJdJD8qVhB9h8BELRQI3VGDZqDSG7Pq6YidvA3wdA2HhfkjQgWzTyC6W8ak0s0jYiZPlm4UvFt4daSWkog0/+hLAvG0nOB0xCzTG2urRVpGVRLKJvkOU1IEy8r0i8NTSllleB7Sfn7yO35o3JmrSYkLvaVpXxNbWLp9wTCDEWKtVMFu79xjRWPeUBee9icLxf6hplpVtDU7jyqtDaqMThKx2kLpRptetr4cuyJeygo+DAyBveEvKhD4LifSkfEtXMvjVMYHlV1p2j4udSXoljGhUlCu5cH8bcV0OIAqQSSGMOs6mO8qWTBuJ9IuVDCJmA8qqqliZimXRC+/VTH6+5P3O5UwinWK1IUQCzSESTT4EiEO/nwGmgnL+zy6tu7Zw4N5pRzy2Kzz1J7DNSNiUxi7ByFKBdJKLLt1RDO94vF7KP8BaIkJK/YwmV5DP6ZzL3OKlak3UrfyM0dEtRgJdz1l4Dkl1LM97PTm8ksUJsAxBS8ncTo7za/v2Y2pC0aMyar1OZNvL8KRAF0OQtIObILWpFeuO0P192AvmQmr8DrcHvJ4TGbTv22oR+j2WlzJ8yowCqvATEHrklrUhPnR0dviiXMVNA4Pzdhnw/OTR2Hi1L0xqamsScQp0/hefqPeQlIGDklkpRK0KukBtbd7ENgfy9ONqBej9XaOyYP1lz6BaLaOD8qbYCu4BKLkUJXgabMvN3KQtBCA6NS2qkIho0f6otLCBgkYgspRS1Il3zaQK1NMjK32U+0MeWizJCofGRD2qk+jMwf6otJCBokYgstRSl3XG4S4Pa+fvFMSKflRmYCwNDY3HgEoto7PlTbSEBQYtEZIGlqOIyd2lQO3+vjlorzQcA7ybn72BoLBPKwW23/qCQgLCdHFCK+i6sDFMaFAhJfBybGFGQkr8zQmM5XpcIoerPNGEAiVU4uJMTEiInWIoyuggBpcHqqCiRD9NIzd81SSvxug/hoVsIQFKoAXZyUkIEl6JkQmBpUPQyBh9xVgWRvyvxuvfhoVswIFeoAXVyckLE+FeVyHiNYeWM0qA6yFVHrdbmI82qtMH5uxqvex0eusUGVGpGV+GAhMhGKiP2SEcHYpzapUH3IMeYrnXNqiDyd9SkHFtsQGLpCww1iBRAsxMiMTJWCI3VxHO/ssE9/6U9XavMqmwE83fUpBxbgIuJpS+4CicH0IyESI78ZEK3w7W6DWl+BzHI4Va5SkJNyjEF9UHSRDcQargCaO2ESImMKyM2O2sStWbNXfM78CCHn1Vhxus4MQGJqadMiDEKHLwDzeW6I+ObCaERR7XMlPmdC+Agh51VIax4HSkWINlzREKsKtz6uVAA3T4ytmqHfur8DnP+SxRmlSswKYcWA5CSejJLXzWkOe4cFECjIuN28zus+qoscKkVotSGEwOQ6jmMpQDl72XZSz5sgwJoTGSMnN9RBJViEaU9lDQBFTRiVkE4U2ZFH/98OxhAA5GxImh+R35P5qvSUhACl9qw0gIkLrUBV0EcMpLvJ9+YkRIGbznBjIxVAfM7ktiu41oKEjhpAJKX2kCec3HyJzUb1pKccM2RyS3WQka32PM7otiuoy4FCZi07iBpko3tOWL4vDty06RrBCo1Sk7BfAQMFs513EtBAiUqoPxtmKU2UhBQlzUiEeovAzGeIF2n3VKQAIkGqGS+uAgbqKIoQQAp26Nlgp26QQjrOrilIN6IAqjiwxJj5DbmJFvxPcRSSOzUDUYI15HLh16FChh5fvCImZDtWwgxR29jTbJFWRD1Z+TUDUqw66jL9VGhAl4egMT6RlG60KkUhLAC+Q/EAgs7CMBP3WCEcB2lfIgJFbxQB0Al18X6hnPGGkIqF7NS+HgxR2YGAd5M3bCFdR21fAiHCt6oPaC28PDrYm5qTZu5P/nsOo1PNxOrogli9soMApBTN7DwrhOA8iFFHe6gVYkREqG23NVnG5I13EOs1yTP2lYMTAHhElSEMK4jPRoRgPIhRR0AWSfulwkRe3r8Ka1fEQjN2LBk8rD4hcz/FS5BhYVwHdejEf6XDynq2El/s9koEbpH7jKKUaVh4wQXa7qUz24amaCCAl1HeTTC//IhRQqg4gni0tSGiQ0iIegpodJQ1jCBm7pBCec66qMRfpcPKVIAFY1InLi3hWxdT4wR8PMF7tkZ5ou4BFVbSNdRi/jV/pYPKVJdrGhc7rKYAw0xVtQqx9IwhhVm6gYjlOsIoYkXRXzv5e6DisaV3f48NnU17ve0+kvso0sYYVxHCk28KOJ7rXadtECIVC0A9pwBhH90CSGU64ihCa6I75vaj2IiIf/kzaNLbAmeg3MdOTRBFPF9VIdhvmicn0E68tElWJLn4FxHDk0w66l9Usc4CLMHJku4R5cwkoJ6wHU6hCYPS4Fd5Rq4epXsOWzX8S408VEB3kcxcPUql+cwXce70MQ3BXyjSf/rVV54jnehiU8KOCD/61XeeE4gQhO2Ar9Vqf/1Km88x//QBFCn3MvVG8/xOzQB1CkBeeU5/oYmgDonoIfvOWh1UkAP3XPQ6qyAHrbnoNVpAXUWcUCAOCBAHBAgDggQBwSIAwLEAQHigABxQIA4IEAcECAOCBAHBIgDAsQBAeKAAHFAgDggQBwQIA4IEAcEiAMCxAEB4oAAcUCAOCBAHBCg/wcPXYdrDIhSHQAAAABJRU5ErkJggg==" />

<!-- rnb-plot-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


### 3B_Cumulative_frequency


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxucGx0X2N1bXVsYXRpdmUgPC0gY3VtdWxhdGl2ZV9mcmVxdWVuY3lfcGxvdChkZl9ncHQ0KSArXG4gIGxhYnMoeSA9IFwiQ29tYmluZWQgZnJlcXVlbmN5XFxub2YgdG9wIDI1IGRpYWdub3Nlc1wiKVxucGx0X2N1bXVsYXRpdmVcbmBgYCJ9 -->

```r
plt_cumulative <- cumulative_frequency_plot(df_gpt4) +
  labs(y = "Combined frequency\nof top 25 diagnoses")
plt_cumulative

3C_Diversity

plt_div <- df_diversity_ci %>% 
  format_criteria() %>% 
  ggplot(aes(x = criteria, y = shannon))+
  geom_point(size = 1)+
  geom_errorbar(aes(ymin = `2.5%`, ymax = `97.5%`), width = 0.3)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = "", y = "Shannon diversity") 
plt_div

3D_Precision

plt_precision <- manual_box_plot(df_precision_stats)
plt_precision

3_Compiled

plt_fig3 <- plot_grid(
  plot_grid(NULL),
  plot_grid(
    plt_rank,
    NULL,
    plt_cumulative,
    nrow = 1, 
    rel_widths = c(3,0.2,2),
    labels = c("A", "", "B"),
    vjust = 0.2
  ),
  plot_grid(NULL),
  plot_grid(
    plt_div,
    NULL,
    plt_precision,
    nrow = 1,
    rel_widths = c(1,0.1,1),
    labels = c("C", "", "D"),
    vjust = 0.2
  ),
  ncol = 1,
  rel_heights = c(0.05,1,0.05,1)
)
plt_fig3
ggsave(plot=plt_fig3,filename=here("figures/4_Diagnosis_diversity.pdf"), width = 6.5, height = 6)

Table

diagnosis_rank_table(df_gpt4, "mast |mastoc|anaphylaxis") %>% mutate(diagnosis = substr(diagnosis, 1, 60)) %>% 
  select(contains(c("diagnosis","mcas"))) %>% 
  head(6) %>% 
  mutate(diagnosis = str_to_sentence(diagnosis)) %>% 
  rename(Diagnosis = diagnosis, `MCAS consortium` = mcas_consortium, `MCAS alternative` = mcas_alternative) %>% 
  flextable() %>% 
  width(width = 2) %>% 
  align(j = 2:3, align = "center", part = "all") %>% 
  {print(., preview = "pdf");.}

Supplemental figures

plot_grid(
  top_diagnosis_plot(df_gpt3, n_diag = 25) + 
    theme(axis.text = element_text(size = 4.5, lineheight = 0.7), 
          strip.text = element_text(size = 7),
          axis.title = element_text(size = 9)) + 
    tidytext::scale_x_reordered(labels = custom_labeler),
  
  top_diagnosis_plot(df_gpt4_icd, n_diag = 25) + 
    theme(axis.text = element_text(size = 4, lineheight = 0.7), 
          strip.text = element_text(size = 7),
          axis.title = element_text(size = 9)) + 
    tidytext::scale_x_reordered(labels = ~custom_labeler(., wrap_width = 45)),
  ncol = 1,
  rel_heights = c(0.9,1),
  labels = c("A","B")
) %>% 
  {ggsave(plot=.,filename=here("figures/S_Top_diagnoses.pdf"), width = 7, height = 10);.}

plot_grid(
  rank_abundance_plot(df_gpt3) +theme(legend.position = c(0.7,0.7)),
  rank_abundance_plot(df_gpt4_icd) +theme(legend.position = c(0.7,0.7)),
  nrow = 1,
  labels = c("A","B")
) %>% 
  {ggsave(plot=.,filename=here("figures/S_Ranked_abundance.pdf"), width = 7.5, height = 3);.}

Alternative combined figure

n_diagnoses_bar <- 15
n_diagnoses_abundance <- 50
n_diagnoses_cumulative <- 50

apply_text_formatting <- list(theme(
  axis.text = element_text(size = 6),
  axis.title = element_text(size = 9),
  legend.text = element_text(size = 6),
  # legend.spacing.y = unit(0.1, 'cm')
  # legend.height = unit(0.1, 'cm'),
  legend.key.height = unit(0.3, 'cm')
  # legend.key.spacing = unit(0.1, 'cm')
  ))
  
  
plot_grid(
  
  ###
  top_diagnosis_plot(df_gpt4, n_diag = n_diagnoses_bar) + 
    theme(axis.text = element_text(size = 5, lineheight = 0.7), 
          strip.text = element_text(size = 7),
          axis.title = element_text(size = 9)) + 
    tidytext::scale_x_reordered(labels = custom_labeler),
  ###
  NULL,
  plot_grid(
    plot_grid(NULL),
    plot_grid(
      NULL,
      ####
      rank_abundance_plot(df_gpt4, n_diagnoses = n_diagnoses_abundance) +theme(legend.position = c(0.7,0.7)) + apply_text_formatting,
      ###
      NULL,
      ###
      cumulative_frequency_plot(df_gpt4, n_diagnoses = n_diagnoses_cumulative) + labs(y = str_glue("Combined frequency\nof top {n_diagnoses_cumulative} diagnoses")) + apply_text_formatting,
      ###
      NULL,
      nrow = 1, 
      rel_widths = c(0.2, 2.5, 0.2, 2, 0.2 ),
      # labels = c("A", "", "B"),
      vjust = 0.2
    ),
    plot_grid(NULL),
    plot_grid(
      NULL,
      ###
       df_diversity_ci %>% 
        format_criteria() %>% 
        ggplot(aes(x = criteria, y = shannon))+
        geom_point(size = 1)+
        geom_errorbar(aes(ymin = `2.5%`, ymax = `97.5%`), width = 0.3)+
        theme_bw() +
        theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
        labs(x = "", y = "Shannon\ndiversity") + apply_text_formatting,
      ###
      NULL,
      ###
      manual_box_plot(df_precision_stats) + labs(y='Bray-Curtis\nsimilarity') + apply_text_formatting,
      ###
      NULL,
      nrow = 1,
      rel_widths = c(0.2, 1, 0.1, 1, 0.2),
      # labels = c("C", "", "D"),
      vjust = 0.2
    ),
    ncol = 1,
    rel_heights = c(0.05,1,0.05,1)
  ),
  
  ncol = 1,
  rel_heights = c(0.9,0.02,1)
) %>% 
  {ggsave(plot=.,filename=here("figures/3_4_combined.pdf"), width = 7, height = 9);.}

n_diagnoses_bar <- 15
n_diagnoses_abundance <- 50
n_diagnoses_cumulative <- 50

apply_text_formatting <- list(theme(
  axis.text = element_text(size = 6),
  axis.title = element_text(size = 9),
  legend.text = element_text(size = 6),
  # legend.spacing.y = unit(0.1, 'cm')
  # legend.height = unit(0.1, 'cm'),
  legend.key.height = unit(0.3, 'cm')
  # legend.key.spacing = unit(0.1, 'cm')
  ))
  
  
plot_grid(
  
  ###
  top_diagnosis_plot(df_gpt4, n_diag = n_diagnoses_bar) + 
    theme(axis.text = element_text(size = 5, lineheight = 0.7), 
          strip.text = element_text(size = 7),
          axis.title = element_text(size = 9)) + 
    tidytext::scale_x_reordered(labels = custom_labeler),
  ###
  NULL,
  plot_grid(
      rank_abundance_plot(df_gpt4, n_diagnoses = n_diagnoses_abundance) +theme(legend.position = c(0.7,0.7)) + apply_text_formatting,
      cumulative_frequency_plot(df_gpt4, n_diagnoses = n_diagnoses_cumulative, width = 0.75) + labs(y = str_glue("Combined frequency\nof top {n_diagnoses_cumulative} diagnoses")) + apply_text_formatting,
      df_diversity_ci %>% 
        format_criteria() %>% 
        ggplot(aes(x = criteria, y = shannon))+
        geom_point(size = 1)+
        geom_errorbar(aes(ymin = `2.5%`, ymax = `97.5%`), width = 0.3)+
        theme_bw() +
        theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
        labs(x = "", y = "Shannon diversity") + apply_text_formatting,
      manual_box_plot(df_precision_stats, width = 0.5) + labs(y='Bray-Curtis similarity') + apply_text_formatting,
      nrow = 1, 
      axis = 'tb',
      align = 'h',
      rel_widths = c(1, 0.7, 0.7, 0.7),
      labels = c(LETTERS[2:5]),
      vjust = 0.2
    ),
  ncol = 1,
  rel_heights = c(1, 0.05, 0.65),
  labels = c("A","","")
)  

Error in count(., criteria, diagnosis, sort = T) : 
  object 'df_gpt4' not found

Alternative final for ICD

n_diagnoses_bar <- 10
n_diagnoses_abundance <- 50
n_diagnoses_cumulative <- 50

title_size <- 9
label_size <- 6
legend_x_pad <- 4
legend_y_pad <- 2

apply_text_formatting <- list(theme(
  axis.text = element_text(size = label_size),
  axis.title = element_text(size = title_size),
  legend.text = element_text(size = label_size),
  strip.text = element_text(size = label_size+1),
  legend.key.height = unit(0.4, 'cm'),
  legend.box.background = element_rect(color = "black", size = 1),
  legend.margin = margin(t = legend_y_pad, r = legend_x_pad, b = legend_y_pad, l = legend_x_pad*1.1),
  legend.spacing.x = unit(0, 'cm'),                           # Horizontal spacing between legend items
  # legend.spacing.y = unit(0, 'cm'),
  # legend.box.spacing = unit(0, "cm")
  ))

strip_margin <- 1
strip_formatting <- list(theme(
  strip.text.x = element_text(margin = margin(t = strip_margin, r = strip_margin, b = strip_margin, l = strip_margin)), 
  strip.text.y = element_text(margin = margin(t = strip_margin, r = strip_margin, b = strip_margin, l = strip_margin))
  # strip.background = element_rect(margin = margin(t = strip_margin, r = strip_margin, b = strip_margin, l = strip_margin))
))

plt_diags <- multi_top_diagnosis_plot(distribution_vis = "points", wrap_width=58, n_diag = n_diagnoses_bar,
                         df_gpt3.5_icd, df_gpt4.0_icd, df_claude3_haiku_t1.0_icd, 
                         df_claude3_opus_t1.0_icd, df_gemini1.0_pro_t1.0_icd) +
  theme(legend.position = "bottom", legend.direction = "horizontal") +
  apply_text_formatting +
  theme(axis.text.y = element_text(size = 6.5)) +
  strip_formatting +
  # theme(legend.position = c(-1,0))+
  theme(panel.spacing = unit(0, "lines")) +
  guides(
  color = guide_legend(override.aes = list(size = 2))  # Increase the point size in the legend
)
  

plt_rank <- multi_ranked_abundance_plot(df_gpt3.5_icd, df_gpt4.0_icd, df_claude3_haiku_t1.0_icd, 
                            df_claude3_opus_t1.0_icd, df_gemini1.0_pro_t1.0_icd)+
  theme(legend.position = "bottom", legend.direction = "horizontal") +
  apply_text_formatting+
  guides(color = guide_legend(ncol = 2))

plt_cumulative <- multi_cumulative_frequency_plot(n_diagnoses = n_diagnoses_cumulative, distribution_vis = "points",
                                                  df_gpt3.5_icd, df_gpt4.0_icd, df_claude3_haiku_t1.0_icd, 
                                df_claude3_opus_t1.0_icd, df_gemini1.0_pro_t1.0_icd) +
  theme(legend.position = "bottom", legend.direction = "horizontal") +
  apply_text_formatting+
  guides(color = guide_legend(ncol = 2)) +
  labs(y = "Combined frequency\nof top 50 diagnoses", x = NULL)

plt_shannon <- tibble(model = c("ChatGPT 3.5", "ChatGPT 4.0", "Claude3 Haiku", "Claude3 Opus", "Gemini 1.0 Pro"),
       data = list(df_gpt3.5_icd, df_gpt4.0_icd, df_claude3_haiku_t1.0_icd, df_claude3_opus_t1.0_icd, df_gemini1.0_pro_t1.0_icd)) %>% 
  mutate(shannon = map(data, calculate_shannon)) %>% 
  select(model, shannon) %>% 
  unnest_wider(shannon) %>% 
  pivot_longer(-model, names_to = "criteria", values_to = "shannon") %>% 
  format_criteria() %>% 
  ggplot(aes(x = criteria, y = shannon))+
  stat_summary(fun.y = mean, aes(group = criteria), geom = "crossbar", size = 0.3, width = 0.75) +
  geom_point(aes(color = model), size = 0.5, position = position_dodge(width = 0.75))+
  theme_bw() +
  scale_color_brewer(palette = "Dark2")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  ggpubr::geom_pwc(aes(group = criteria), 
                   method = "wilcox.test", 
                   p.adjust.method = "BH", 
                   hide.ns = T, 
                   label = "p.adj.signif", 
                   bracket.nudge.y = 0.1, 
                   vjust = 0.6, 
                   step.increase = 0.14, 
                   tip.length = 0.02
                   )+
  labs(x = NULL, y = "Shannon Diversity", color = "")+
  apply_text_formatting+
  theme(legend.position = "bottom", legend.direction = "horizontal") +
  guides(color = guide_legend(ncol = 2))

plt_precision <- df_precision %>% 
  format_criteria() %>% 
  format_models() %>%
  ggplot(aes(x = criteria, y = y))+
    stat_summary(fun.y = mean, aes(group = criteria), geom = "crossbar", size = 0.3, width = 0.75) +
  geom_point(aes(color = model), size = 0.5, position = position_dodge(width = 0.75))+
  theme_bw()+
  theme(axis.text.x = element_text(angle= 45, hjust = 1))+
  labs(x=NULL, y = "Average Bray-Curtis\nSimilarity") +
  ggpubr::geom_pwc(method = "wilcox.test", p.adjust.method = "BH", hide.ns = T, label = "p.adj.signif", bracket.nudge.y = 0.3, vjust = 0.6, step.increase = 0.14, tip.length = 0.02) +
  ylim(c(0,1)) +
  labs(color=NULL)+
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "none")+
  apply_text_formatting




full_plt <- plot_grid(
  
  ###
  plt_diags,
  ###
  NULL,
  plot_grid(
      plt_rank,
      plt_cumulative,
      plt_shannon,
      plt_precision,
      nrow = 1, 
      axis = 'tb',
      align = 'h',
      rel_widths = c(1, 0.7, 0.7, 0.7),
      labels = c(LETTERS[2:5]),
      vjust = 0.2
    ),
  ncol = 1,
  rel_heights = c(1.2, 0.05, 0.65),
  labels = c("A","","")
)  

full_plt

Things to fix - Legend position for C-E - Legend width for B - Move legend for A to the left of “Frequency?” - Rank plot line weight

ggsave(plot=full_plt,filename=here("figures/3_diagnosis_diversity.pdf"), width = 7.5, height = 7.5)

set_table_properties(opts_pdf = list(tabcolsep = 0))

set_flextable_defaults(fonts_ignore=TRUE)

multi_diagnosis_rank_table(search_pattern = "T78\\.2 |D47\\.02 |D89\\.41 |D89\\.49 |D89\\.4 ",
                                         df_gpt3.5_icd, df_gpt4.0_icd, df_claude3_haiku_t1.0_icd, df_claude3_opus_t1.0_icd, df_gemini1.0_pro_t1.0_icd) %>% 
  flextable() %>% 
  width(width = 2) %>% 
  fontsize(size = 9) %>% 
  fontsize(size = 10, part = "header") %>% 
  padding(padding = 0) %>% 
  align(j = 2:3, align = "center", part = "all") %>% 
  set_table_properties(opts_pdf = list(arraystretch = 1.25)) %>% 
  {print(., preview = "pdf");.}

Diagnosis	MCAS - Consortium	MCAS - Alternative
T78.2 Anaphylactic shock, unspecified	1 [1, 1, 1, 1, 1]	151 [216, 87, 96, 186, 168]
D47.02 Systemic mastocytosis	8 [22, 8, 6, 2, 2]	56 [92, 76, 24, 45, 41]
D89.41 Monoclonal mast cell activation syndrome	40 [128, 23, 28, 11, 11]	61 [141, 68, 23, 39, 36]
D89.49 Other mast cell activation disorder	152 [307, 63, 96, 137, 155]	661 [1156, 557, 176, 383, 1031]
D89.4 Mast cell activation syndrome and related disorders	452 [725, 178, NA, NA, NA]	1380 [NA, 870, 1845, 1424, NA]